Contents 



On 
On 

c 
3 



(N 



> 

m 
en 
m 

N£> 
O 
ON 
On 

0,; 

6 



13 



1 INTRODUCTION 



O Current fields of research 



L2 Overview of the numerical methods 



1.3 Plan of the Review 



12 SPECIAL RELATIVISTIC HYDRO-DYNAMICSI 



2.1 Equations 



2.2 SRHD as a hyperbolic system of conservation laws 



2.3 Exact solution of the Riemann problem in SRHD 



3 HIGH-RESOLUTION SHOCK-CAPTURING 



METHODS 



3 

3 
4 

5 

5 

5 

7 
7 



11 



3.1 Rclativistic PPM| 11 

3~2 Relativistic Glimm's method! 12 



3.3 Two-shock approximation for relativistic hydro dynamics| 12 



3.4 Roe-type relativistic solvers| 14 

3.5 Falle and Komissarov upwind scheme] 15 

3~6 Relativistic HLL Method! 16 



3.7 Marquina's flux formula 17 

Symmetric TVD schemes with nonlinear numerical dissipation! 18 



3.8 



4 OTHER DEVELOPMENTS 



18 

4.1 Van Putten's approach! 18 

19 

22 



4.2 Relativistic SPH 



4~3" Relativistic beam scheme 



E SUMMARY OF METHODS! 



16 TEST BENCH! 



23 



23 



3.1 Relativistic shock heating in planar, cylindrical and spherical geometry! . . 23 

3.2 Propagation of relativistic blast waves| 28 

32 



p. 2.1 ProblemT 



|3X2~ 



Problem 2| 35 

40 



p.2.3 Collision of two relativistic blast waves 



7 APPLICATIONS 



44 

7.1 Astrophysical jets| 44 

7.2 Gamma-Ray Bursts (GRBs]| 49 



8 CONCLUSION 



51 

O Evaluation of the methods' 51 



3.2 Further developments! 52 



p. 2.1 incorporation of realistic microphysicsj 52 

P.2.2 Coupling of SRHD schemes with AMR] 54 



p. 2. 3 General rclativistic hydrodynamics (GliHD) 54 

p. 2.4 Relativistic magneto-hydrodynamics (RMHD)| 55 



9 ADDITIONAL INFORMATION 



3.1 Algorithms to recover primitive quantities 



56 

56 

3.2 Spectral decomposition of the 3D SRHD equations! 57 

-J. 3 Program RIEMANN] 58 

3.4 Basics of HRSC methods and recent developments! 59 

3.5 Newtonian SPH equations! 61 



NUMERICAL HYDRODYNAMICS IN SPECIAL 

RELATIVITY 

J. M. Marti 

Departamento de Astronomia y Astrofisica 

Universidad de Valencia 

46100 Burjassot (Valencia), Spain 

E. Miiller 

Max-Planck-Institut fiir Astrophysik 
Karl-Schwarzschild-Str. 1, 85748 Garching, Germany 

February 1, 2008 

ABSTRACT 

This review is concerned with a discussion of numerical methods for the solution of the 
equations of special relativistic hydrodynamics (SRHD). Particular emphasis is put on 
a comprehensive review of the application of high-resolution shock-capturing methods 
in SRHD. Results obtained with different numerical SRHD methods are compared, and 
two astrophysical applications of SRHD flows are discussed. An evaluation of the various 
numerical methods is given and future developments are analyzed. 

1 INTRODUCTION 

1.1 Current fields of research 

Relativity is a necessary ingredient for describing astrophysical phenomena involving com- 
pact objects. Among these phenomena are core collapse supernovae, X-ray binaries, pul- 
sars, coalescing neutron stars, formation of black holes, micro-quasars, active galactic 
nuclei, superluminal jets and gamma-ray bursts. General relativistic effects must be con- 
sidered when strong gravitational fields are encountered as, for example, in the case of 
coalescing neutron stars or near black holes. The significant gravitational wave signal pro- 
duced by some of these phenomena can also only be understood in the framework of general 
theory of relativity. There are, however, astrophysical phenomena which involve flows at 



relativistic speeds but no strong gravitational fields and thus at least certain aspects of 
these phenomena can be described within the framework of special relativity. 

Another field of research, where special relativistic "flows" are encountered, are present- 
day heavy-ion collision experiments taking place in large particle accelerators. The heavy 
ions are accelerated to ultra-relativistic velocities very close to the speed of light (~ 99.998% 
162| ) to study the equation of state for hot dense nuclear matter. 



1.2 Overview of the numerical methods 

The first attempt to solve the equations of relativistic hydrodynamics (SRHD) was made 
by Wilson 1 183 , 184 ] and collaborators [27, 75 1 using an Eulerian explicit finite difference 



code with monotonic transport. The code relies on artificial viscosity techniques ||181 
15 Of to handle shock waves. It has been widely used to simulate flows encountered in 
cosmology, axisymmetric relativistic stellar collapse, accretion onto compact objects and, 
more recently, collisions of heavy ions. Almost all the codes for numerical (both special and 
general) relativistic hydrodynamics developed in the eighties ||139| , |163| , |123| , |122| , |124| , |52| 
were based on Wilson's procedure. However, despite its popularity it turned out to be 
unable to accurately describe extremely relativistic flows (Lorentz factors larger than 2; 
see, e.g. 



In the mid eighties, Norman & Winkler |128|1 proposed a reformulation of the difference 
equations with artificial viscosity consistent with the relativistic dynamics of non-perfect 
fluids. The strong coupling introduced in the equations by the presence of the viscous terms 
in the definition of relativistic momentum and total energy densities required an implicit 
treatment of the difference equations. Accurate results across strong relativistic shocks 
with large Lorentz factors were obtained in combination with adaptive mesh techniques. 
However, no multidimensional version of this code was developed. 

Attempts to integrate the SRHD equations avoiding the use of artificial viscosity were 
performed in the early nineties. Dubai [-1 (] developed a 2D code for relativistic magneto- 
hydrodynamics based on an explicit second-order Lax-Wendroff scheme incorporating a 
flux-corrected transport (FCT) algorithm [[jl|. Following a completely different approach 
Mann [|9JJ proposed a multidimensional code for general relativistic hydrodynamics based 
on smoothed particle hydrodynamics (SPH) techniques |119|| , which he applied to rela- 
tivistic spherical collapse |100|| . When tested against ID relativistic shock tubes all these 
codes performed similar to the code of Wilson. More recently, Dean et al. |38| have applied 
flux correcting algorithms for the SRHD equations in the context of heavy ion collisions. 

L2I1 



Recent developments in relativistic SPH methods |29|, |160| are discussed in Section 

A major break-through in the simulation of ultra-relativistic flows was accomplished 
when high-resolution shock-capturing (HRSC) methods, specially designed to solve hyper- 
bolic systems of conservations laws, were applied to solve the SRHD equations ||104| , |103| , 
50, 5jJ. This review is intended to provide a comprehensive discussion of different HRSC 



methods and of related methods used in SRHD. Numerical methods for special relativistic 
MHD flows are not included, because they are beyond the scope of this review. However, 
we may include such a discussion in a future update of this article. 



1.3 Plan of the Review 

The review is organized as follows. Section (|2]) contains a derivation of the equations of 
special relativistic (perfect) fluid dynamics, as well as a discussion of their main proper- 
ties. In Section (|3]) the most recent developments in numerical methods for SRHD are 
reviewed paying particular attention to high-resolution shock-capturing methods. Other 
developments in special relativistic numerical hydrodynamics are discussed in Section (^). 
Numerical results obtained with different methods as well as analytical solutions for sev- 
eral test problems are presented in Section ([]). Two astrophysical applications of SRHD 
are discussed in Section (|7|). An evaluation of the various numerical methods is given in 
Section (^) together with an outlook for future developments. Finally, some additional 
technical information is presented in Section (pf). 



The reader is assumed to have basic knowledge in classical |91| , p4J1 and relativistic fluid 



dynamics |166| , |6|, as well as in finite difference/volume methods for partial differential 



equations ||148 , |1 29| . A discussion of modern finite volume methods for hyperbolic systems 



of conservation laws can be found, e.g., in [95, 96, B2J. The theory of spectral methods 



for fluid dynamics is developed in | 23fl , and smoothed particle hydrodynamics (SPH) is 



reviewed in |L19|. 



2 SPECIAL RELATIVISTIC HYDRO- DYNAMICS 

2.1 Equations 

Using the Einstein summation convention the equations describing the motion of a rela- 
tivistic fluid are given by the five conservation laws 

(P« M );m = 0, (1) 

I*£ = , (2) 

where (p, v = 0, ...,3), and where ; p denotes the covariant derivative with respect to 
coordinate x M . Furthermore, p is the proper rest-mass density of the fluid, u^ its four- 
velocity, and T^ v is the stress-energy tensor, which for a perfect fluid can be written as 

T" u = phuV + pg^ . (3) 

Here, g^ u is the metric tensor, p the fluid pressure, and and h the specific enthalpy of the 
fluid defined by 

h = 1 + e + p/p , (4) 

where e is the specific internal energy. Note that we use natural units (i.e., the speed of 
light c = 1) throughout this review. 

In Minkowski space-time and Cartesian coordinates (t,x x ,x 2 ,x 3 ), the conservation 
equations ([!]), (|2|) can be written in vector form as 

du <9F*(u) 

^ + ^^ = ' (5) 



where % = 1,2, 3. The state vector u is defined by 

u=(D,S\S 2 ,S 3 ,t) t (6) 

and the flux vectors F* are given by 

F l = {Dv\ SV + P 5 U , SV + p5 2 \ SV + P 5 3 \ S l - Drff . (7) 

The five conserved quantities D, S l ,S 2 ,S 3 and r are the rest-mass density, the three 
components of the momentum density, and the energy density (measured relative to the 
rest mass energy density), respectively. They are all measured in the laboratory frame, and 
are related to quantities in the local rest frame of the fluid (primitive variables) through 

D = pW, (8) 

& = phW 2 v l i = 1, 2, 3 , (9) 

r = phW 2 - p - D , (10) 

where v % are the components of the three-velocity of the fluid 

v l = u l /u° (11) 

and W is the Lorentz factor 

W = u°= l % . (12) 

yl — v % Vi 

The system of equations @ with definitions @ - (0) is closed by means of an equation 
of state (EOS), which we shall assume to be given in the form 

p = p(p,e). (13) 

In the non-relativistic limit (i.e.,v <C 1, h — > 1) D, S l and r approach their Newtonian 
counterparts p, pv % and pE = pe + pv 2 /2, and equations (H) reduce to the classical ones. In 
the relativistic case the equations of system (|) are strongly coupled via the Lorentz factor 



and the specific enthalpy, which gives rise to numerical complications (see Section |2T3| ). 

In classical numerical hydrodynamics it is very easy to obtain v l from the conserved 
quantities (i.e.,p and pv l ). In the relativistic case, however, the task to recover (p,v\p) 
from (D, S l , t) is much more difficult. Moreover, as state-of-the-art SRHD codes are based 
on conservative schemes where the conserved quantities are advanced in time, it is necessary 
to compute the primitive variables from the conserved ones one (or even several) times per 
numerical cell and time step making this procedure a crucial ingredient of any algorithm. 



(see Sect. |97T 



2.2 SRHD as a hyperbolic system of conservation laws 

An important property of system ([5]) is that it is hyperbolic for causal EOS 0. For 
hyperbolic systems of conservation laws, the Jacobians <9F*(u)/<9u have real eigenvalues and 



a complete set of eigenvectors (see Sect. |9.2|) . Information about the solution propagates at 
finite velocities given by the eigenvalues of the Jacobians. Hence, if the solution is known 
(in some spatial domain) at some given time, this fact can be used to advance the solution 
to some later time (initial value problem). However, in general, it is not possible to derive 
the exact solution for this problem. Instead one has to rely on numerical methods which 
provide an approximation to the solution. Moreover, these numerical methods must be able 
to handle discontinuous solutions, which are inherent to non-linear hyperbolic systems. 

The simplest initial value problem with discontinous data is called a Riemann problem, 
where the one dimensional initial state consists of two constant states separated by a 
discontinuity. The majority of modern numerical methods, the so-called Godunov-type 
methods, are based on exact or approximate solutions of Riemann problems. Because of 
its theoretical and numerical importance, we discuss the solution of the special relativistic 
Riemann problem in the next subsection. 

2.3 Exact solution of the Riemann problem in SRHD 

Let us first consider the one dimensional special relativistic flow of an ideal gas with an 
adiabatic exponent 7 in the absence of a gravitational field. The Riemann problem then 
consists of computing the breakup of a discontinuity, which initially separates two arbitrary 
constant states L (left) and R (right) in the gas (see Fig. [I] with L = 1 and R = 5). For 
classical hydrodynamics the solution can be found, e.g., in fl34|. In the case of SRHD, the 



Riemann problem has been considered by Marti & Miiller ||105| , who derived an exact 
solution generalising previous results for particular initial data |168|| . 

The solution to this problem is self-similar, because it only depends on the two constant 
states defining the discontinuity vl and vr, where v = (p,p,v), and on the ratio (x — 
xo)/(t — to), where xq and to are the initial location of the discontinuity and the time of 
breakup, respectively. Both in relativistic and classical hydrodynamics the discontinuity 
decays into two elementary nonlinear waves (shocks or rarefactions) which move in opposite 
directions towards the initial left and right states. Between these waves two new constant 
states v L;(i and v R * (note that v L * = 3 and vr* = 4 in Fig. [[]) appear, which are separated 
from each other through a contact discontinuity moving with the fluid. Across the contact 
discontinuity the density exhibits a jump, whereas pressure and velocity are continuous (see 
Figure [l]). As in the classical case, the self-similar character of the flow through rarefaction 
waves and the Rankine-Hugoniot conditions accross shocks provide the relations to link 
the intermediate states vs* (S = L, R) with the corresponding initial states vg. They also 
allow one to express the fluid flow velocity in the intermediate states v$* as a function of 
the pressure ps* in these states. Finally, the steadiness of pressure and velocity across the 
contact discontinuity implies 

v L *(p*) =fR*(P*) 5 (14) 



where p* = pl* = Pr*, which closes the system. The functions vs*(p) are defined by 



where TZ s (p) (S s (p)) denotes the family of all states which can be connected through a 
rarefaction (shock) with a given state vs ahead of the wave. 

The fact that one Riemann invariant is constant through any rarefaction wave provides 
the relation needed to derive the function TZ S 



TZ s (p) 



'l + v s )A±(p) 



vs 



(l + v 3 )A±(p) + (l-v s ) 



with 



A±(p) 



V7- 1 - c (p) V7- 1 + cs 



±- 



(16) 



(17) 



v V7 - 1 + c(p) V7- 1 ~ c s, 
the + (— ) sign of A± corresponding to S — L (S — R). In the above equation, cs is the 
sound speed of the state V5, and c(p) is given by 



\(7- l)ps(j>/Ps) lh + 1P) 

The family of all states S s (p), which can be connected through a shock with a given 
state vs ahead of the wave, is determined by the shock jump conditions. One obtains 



S s (p) 



h s W s v s ± 



P-Ps 



j(p)yJl-V±(p)\ 



h s W s + (p- ps) 



± 



v s 



psW s j( P )Ji-v±( P y 



-1 



(19) 



where the + (— ) sign corresponds to S = R (S = L). V±(p) and j(p) denote the shock 
velocity and the modulus of the mass flux across the shock front, respectively. They are 
given by 



V±(p) 



p s wiv s ±j(p) 2 ^i + (ps/j(p)y 



pWs+j(py 



and 



jip) 



\ 



Ps-P 



hi 



h(pf 



2h.c 



(20) 
(21) 



Ps - P Ps 

where the enthalpy h(p) of the state behind the shock is the (unique) positive root of the 
quadratic equation 



(7-l)(Ps-p) 



IP 



h2 _(l-l)(ps- P ) h+ hs(ps-p)_ h2 







IP 

8 



Ps 



(22) 
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Figure 1: Schematic solution of a Riemann problem in special relativistic hydrodynamics. The 
initial state at t = (top figure) consists of two constant states (1) and (5) with p\ > p^, p\ > ps 
and V\ = t>2 = separated by a diaphragm at x-q. The evolution of the flow pattern once the 
diaphragm is removed (middle figure) is illustrated in a space-time diagram (bottom figure) with 
a shock wave (solid line) and a contact discontinuity (dashed line) moving to the right. The 
bundle of solid lines represents a rarefaction wave propagating to the left. 
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Figure 2: Graphical solution in the p-v plane of the Riemann problems defined by the initial 
states (j?l = 10 3 , pl = 1, ut, = 0.5) and (p l R , /9r = 1, ur = 0) with i = 1,2,3,4, where p^ = 10 2 , 
p^ = 10, p^ = 1, and p^ = 10 _1 , respectively. The adiabatic index of the fluid is 5/3 in all cases. 
Note the asymptotic behavior of the functions as they approach v = 1 (i.e., the speed of light). 



which is obtained from the Taub adiabat (the relativistic version of the Hugoniot adiabat) 
for an ideal gas equation of state. 

The functions vl*(p) and ur*(p) are displayed in Fig. |] in a p-v diagram for a particular 
set of Riemann problems. Once p* has been obtained, the remaining state quantities and 
the complete Riemann solution, 



u = u(- — — ;u L ,u R )), 

t — In 



(23) 



can easily be derived. 



In Section p73| we provide a FORTRAN program called RIEMANN, which allows one to 
compute the exact solution of an arbitrary special relativistic Riemann problem using the 
algorithm just described. 

The treatment of multidimensional special relativistic flows is significantly more difficult 
than that of multidimensional Newtonian flows. In SRHD all components (normal and 
tangential) of the flow velocity are strongly coupled through the Lorentz factor, which 
complicates the solution of the Riemann problem severely. For shock waves, this coupling 
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'only' increases the number of algebraic jump conditions, which must be solved. However, 



for rarefactions it implies the solution of a system of ordinary differential equations [105] - 



3 HIGH-RESOLUTION SHOCK-CAPTURING 
METHODS 

The application of high-resolution shock-capturing (HRSC) methods caused a revolution 
in numerical SRHD. These methods satisfy in a quite natural way the basic properties 
required for any acceptable numerical method: (i) high order of accuracy, (ii) stable and 
sharp description of discontinuities, and (iii) convergence to the physically correct solu- 
tion. Moreover, HRSC methods are conservative, and because of their shock capturing 
property discontinuous solutions are treated both consistently and automatically when- 
ever and wherever they appear in the flow. 

As HRSC methods are written in conservation form, the time evolution of zone aver- 
aged state vectors is governed by some functions (the numerical fluxes) evaluated at zone 
interfaces. Numerical fluxes are mostly obtained by means of an exact or approximate 
Riemann solver. High resolution is usually achieved by using monotonic polynomials in 
order to interpolate the approximate solutions within numerical cells. 

Solving Riemann problems exactly involves time-consuming computations, which are 
particularly costly in the case of multidimensional SRHD due to the coupling of the equa- 
tions through the Lorentz factor (see Section p3|) . Therefore, as an alternative, the usage 
of approximate Riemann solvers has been proposed. 

In this Section we summarize the computation of the numerical fluxes in a number of 
methods for numerical SRHD. Methods based on exact Riemann solvers are discussed in 
Sections[3J] and |3]^, while those based on approximate solvers are discussed in Sections ( p73 



Readers not familiar with HRSC methods are referred to Sectionl^J, where the basic 
properties of these methods as well as an outline of the recent developments are described. 

3.1 Relativistic PPM 



Marti & Miiller ||106| have used the procedure discussed in Section|2]3| to construct an exact 



Riemann solver, which they then incorporated in an extension of the PPM method |H| for 
ID SRHD. In their relativistic PPM method numerical fluxes are calculated according to 

F RPPM = F(u(0;u L ,u R )), (24) 

where u^ and u# are approximations of the state vector at the left and right side of a 
zone interface obtained by a second-order accurate interpolation in space and time, and 
u(0; ul, ur) is the solution of the Riemann problem defined by the two interpolated states 
at the position of the initial discontinuity. 

The PPM interpolation algorithm described in p2] gives monotonic conservative parabolic 



profiles of variables within a numerical zone. In the relativistic version of PPM the orig- 
inal interpolation algorithm is applied to zone averaged values of the primitive variables 

11 



v = (p,p,v), which are obtained from zone averaged values of the conserved quantities u. 
For each zone j, the quartic polynomial with zone-averaged values a_j_2, %-i, cij, %+i, and 
a j+2 (where a = p,p,v) is used to interpolate the structure inside the zone. In particular, 
the values of a at the left and right interface of the zone, a^j and cirj, are obtained this 
way. These reconstructed values are then modified such that the parabolic profile, which 
is uniquely determined by a L j, a^j and a,-, is monotonic inside the zone. 

The time-averaged fluxes at an interface j + 1/2 separating zones j and j + 1 are 
computed from two spatially averaged states, v ■ , 1 T and v.- , 1 R at the left and right 
side of the interface, respectively. These left and right states are constructed taking into 
account the characteristic information reaching the interface from both sides during the 
time step. In the relativistic version of PPM the same procedure as in [^] has been followed 
using the characteristic speeds and Riemann invariants of the equations of relativistic 
hydrodynamics. 

3.2 Relativistic Glimm's method 



Wen et aJ. [ p.82|| have extended Glimm's random choice method [65] to ID SRHD. They 



developed a first-order accurate hydrodynamic code combining Glimm's method (using an 
exact Riemann solver) with standard finite difference schemes. 

In the random choice method, given two adjacent states, u™ and u™ +1 , at time t n , the 
value of the numerical solution at time t n+1 ' 2 and position Xj+ 1/2 is given by the exact 
solution u(x, t) of the Riemann problem evaluated at a randomly chosen point inside zone 
(j,3 + 1), i-e., 

u i+a- u \> + DAt' Ui '^ + v ' [ } 

where £ n is a random number in the interval [0, 1] . 

Besides being conservative on average, the main advantages of Glimm's method are 
that it produces both completely sharp shocks and contact discontinuities, and that it is 
free of diffusion and dispersion errors. 

Chorin [^SJ applied Glimm's method to the numerical solution of homogeneous hyper- 



bolic conservation laws. Colella |30| proposed an accurate procedure of randomly sampling 



the solution of local Riemann problems and investigated the extension of Glimm's method 
to two dimensions using operator splitting methods. 

3.3 Two-shock approximation for relativistic hydrodynamics 

This approximate Riemann solver is obtained from a relativistic extension of Colella's 
method |3(J for classical fluid dynamics, where it has been shown to handle shocks of 



arbitrary strength [RD, |186|. In order to construct Riemann solutions in the two-shock ap- 



proximation one analytically continues shock waves towards the rarefaction side (if present) 
of the zone interface instead of using an actual rarefaction wave solution. Thereby one gets 
rid of the coupling of the normal and tangential components of the flow velocity (see 
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Section p73|) , and the remaining minor algebraic complications are the Rankine-Hugoniot 
conditions across oblique shocks. Balsara |§ has developed an approximate relativistic 
Riemann solver of this kind by solving the jump conditions in the shocks' rest frames 
in the absence of transverse velocities, after appropriate Lorentz transformations. Dai & 
Woodward |35[] have developed a similar Riemann solver based on the jump conditions 



across oblique shocks making the solver more efficient. 

Table |l| gives the converged solution for the intermediate states obtained with both 
Balsara's and Dai's & Woodward's procedure for the case of the Riemann problems defined 
in Section |6^ (involving strong rarefaction waves) together with the exact solution. Despite 
the fact that both approximate methods involve very different algebraic expressions, their 
results differ by less than 2%. However, the discrepancies are much larger when compared 
with the exact solution (up to a 100% error in the density of the left intermediate state 
in Problem 2). The accuracy of the two-shock approximation should be tested in the 
ultra-relativistic limit, where the approximation can produce large errors in the Lorentz 
factor (in the case of Riemann problems involving strong rarefaction waves) with important 
implications for the fluid dynamics. Finally, the suitability of the two-shock approximation 
for Riemann problems involving transversal velocities still needs to be tested. 

Table 1: Pressure p*, velocity v* and densities /9l* (left), /3r* (right) for the intermediate state 
obtained for the two-shock approximation of Balsara H (B) and of Dai & Woodward |J5| (DW) 
compared to the exact solution (Exact) for the Riemann problems defined in Sect. |6.2|. 



Method p* v* p L * p R * 



Problem 1 










B 


1.440E+00 


7.131E-01 


2.990E+00 


5.069E+00 


DW 


1.440E+00 


7.131E-01 


2.990E+00 


5.066E+00 


Exact 


1.445E+00 


7.137E-01 


2.640E+00 


5.062E+00 


Problem 2 










B 


1.543E+01 


9.600E-01 


7.325E-02 


1.709E+01 


DW 


1.513E+01 


9.608E-01 


7.254E-02 


1.742E+01 


Exact 


1.293E+01 


9.546E-01 


3.835E-02 


1.644E+01 
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3.4 Roe-type relativistic solvers 



Linearized Riemann solvers are based on the exact solution of Riemann problems of a mod- 
ified system of conservation equations obtained by a suitable linearization of the original 
system. This idea was put forward by Roe |151|| , who developed a linearized Riemann 
solver for the equations of ideal (classical) gas dynamics. Eulderink et al. [50, 51 1 have 



extended Roe's Riemann solver to the general relativistic system of equations in arbitrary 
spacetimes. Eulderink uses a local linearization of the Jacobian matrices of the system 
fulfilling the properties demanded by Roe in his original paper. 

Let B = <9F/<9u be the Jacobian matrix associated with one of the fluxes F of the 
original system, and u the vector of unknowns. Then, the locally constant matrix B, 
depending on Ul and % (the left and right state defining the local Riemann problem) 
must have the following four properties: 

1. It constitutes a linear mapping from the vector space u to the vector space F. 

2. As u L -> u R -> u, B(u L , u R ) -> B(u). 

3. For any u L , u R , S(u L , u R )(u R - u L ) = F(u R ) - F(u L ). 

4. The eigenvectors of B are linearly independent. 

Conditions 1 and 2 are necessary if one is to recover smoothly the linearized algorithm 
from the nonlinear version. Condition 3 (supposing 4 is fulfilled) ensures that if a single 
discontinuity is located at the interface, then the solution of the linearized problem is the 
exact solution of the nonlinear Riemann problem. 

Once a matrix, B, satisfying Roe's conditions has been obtained for every numerical 
interface, the numerical fluxes are computed by solving the locally linear system. Roe's 
numerical flux is then given by 



^ROE 



F(u L ) + F(ur) 



]T| A (p)| 5 (p) ? (p) 



(26) 



with 



Q 



(P) 



1« ■ (u R - u L ) , (27) 

where \( p \ f&\ and 1^ are the eigenvalues and the right and left eigenvectors of B, 
respectively (p runs from 1 to the number of equations of the system) . 

Roe's linearization for the relativistic system of equations in a general spacetime can 



be expressed in terms of the average state pDL |5T 



w 



W L + Wr 

k L + k R 



with 



w = (ku°, ku 1 , ku 2 , ku 3 , k—r) 

ph 



(28) 



(29) 
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and 

k 2 = V^ph, (30) 

where g is the determinant of the metric tensor g^. The role played by the density p in 
case of the Cartesian non-relativistic Roe solver as a weight for averaging, is taken over 
in the relativistic variant by k, which apart from geometrical factors tends to p in the 
non-relativistic limit. A Riemann solver for special relativistic flows and the generalization 
of Roe's solver to the Euler equations in arbitrary coordinate systems are easily deduced 
from Eulderink's work. The results obtained in ID test problems for ultra-relativistic flows 
(up to Lorentz factors 625) in the presence of strong discontinuities and large gravitational 
background fields demonstrate the excellent performance of the Eulderink-Roe solver |y] . 
Relaxing condition 3 above, Roe's solver is no longer exact for shocks but still produces 
accurate solutions, and moreover, the remaining conditions are fulfilled by a large number of 
averages. The ID general relativistic hydrodynamic code developed by Romero et al. |153 ] 



uses flux formula (|26|) with an arithmetic average of the primitive variables at both sides 
of the interface. It has successfully passed a long series of tests including the spherical 
version of the relativistic shock reflection (see Section |6~T| ). 

Roe's original idea has been exploited in the so-called local characteristic approach 



(see, e.g., ||193|| ). This approach relies on a local linearization of the system of equations by 
defining at each point a set of characteristic variables, which obey a system of uncoupled 
scalar equations. This approach has proven to be very successful, because it allows for 
the extension to systems of scalar nonlinear methods. Based on the local characteristic 
approach are the methods developed by Marquina et al. [ |103|| and Dolezal & Wong pT| , 
which both use high-order reconstructions of the numerical characteristic fluxes, namely 



PHM |D1 and ENO |I| (see Section K 



3.5 Falle and Komissarov upwind scheme 

Instead of starting from the conservative form of the hydrodynamic equations, one can use 
a primitive-variable formulation in quasi-linear form 

dv <9v 

m +A T^°- (31) 

where v is any set of primitive variables. A local linearization of the above system allows 
one to obtain the solution of the Riemann problem, and from this the numerical fluxes 
needed to advance a conserved version of the equations in time. 

Falle & Komissarov 55j have considered two different algorithms to solve the local 



Riemann problems in SRHD by extending the methods devised in f5"3" |. In a first algorithm, 
the intermediate states of the Riemann problem at both sides of the contact discontinuity, 
v^* and v^*, are obtained by solving the system 

v L * = v L + b L r£ , v Rt = v R + 6 R r£ , (32) 
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where r£ is the right eigenvector of «4(vl) associated with sound waves moving upstream 
and r^ is the right eigenvector of *A(v R ) of sound waves moving downstream. The continu- 
ity of pressure and of the normal component of the velocity across the contact discontinuity 
allows one to obtain the wave strengths b^ and b^ from the above expressions, and hence 
the linear approximation to the intermediate state v*(vl, vr). 



In the second algorithm proposed by Falle & Komissarov [55], a linearization of system 
(J31~l) is obtained by constructing a constant matrix -4.(v L ,v R ) = ^4(^(vl + v R )). The 
solution of the corresponding Riemann problem is that of a linear system with matrix A, 
i.e., 

v* = v L +^ 5^?W, (33) 

A(p)<0 



or, equivalent ly, 



with 



v, = v R -_X) 5W?W), (34) 

A(p)>0 



« (ri = l (p) ■ (v R - v L ) , (35) 

where A^, r^ p \ and 1^ are the eigenvalues and the right and left eigenvectors of A, 
respectively (p runs from 1 to the number of equations of the system) . 

In both algorithms, the final step involves the computation of the numerical fluxes for 
the conservation equations 

F FK = F(u(v*(v L ,v R ))). (36) 

3.6 Relativistic HLL Method 

Schneider et al. ||157|| have proposed to use the method of Harten, Lax & van Leer (HLL 




hereafter [74|) to integrate the equations of SRHD. This method avoids the explicit cal- 
culation of the eigenvalues and eigenvectors of the Jacobian matrices and is based on an 
approximate solution of the original Riemann problems with a single intermediate state 

for x < Oht 
u HLL (x/t; u L , Ur) = | u» for a^t < x < a R £ , (37) 

for x > aRt 

where a^ and a R are lower and upper bounds for the smallest and largest signal velocities, 
respectively. The intermediate state u* is determined by requiring consistency of the 
approximate Riemann solution with the integral form of the conservation laws in a grid 
zone. The resulting integral average of the Riemann solution between the slowest and 
fastest signals at some time is given by 

a R u R - a L u L - F(u R ) + F(u L ) 

U* — , [OQ) 

a R - a L 
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and the numerical flux by 



,i ill «r f ( u l) - a h F(u R ) + a£a L (u R - u L ) 



TpHLL _ "'H'- V~W "X^V~n/ ' "K^LV^ ^W /on") 



OjJ> CLj 



where 

a^ = min{0, a L } , a^ = max{0, a R } . (40) 

An essential ingredient of the HLL scheme are good estimates for the smallest and 
largest signal velocities. In the non-relativistic case, Einfeldt [49| proposed calculating 



them based on the smallest and largest eigenvalues of Roe's matrix. This HLL scheme with 
Einfeldt 's recipe is a very robust upwind scheme for the Euler equations and possesses the 
property of being positively conservative. The method is exact for single shocks, but it is 
very dissipative, especially at contact discontinuities. 

Schneider et al. [ |157|| have presented results in ID ultra-relativistic hydrodynamics using 
a version of the HLL method with signal velocities given by 

a R = (v + c 8 )/(l + vc 8 ), (41) 

a h = (v-c t )/(l-vc,), (42) 

where c s is the relativistic sound speed, and where the bar denotes the arithmetic mean 



between the initial left and right states. Duncan & Hughes JJBJ have generalized this 
method to 2D SRHD and applied it to the simulation of relativistic extragalactic jets. 

3.7 Marquina's flux formula 

Godunov-type schemes are indeed very robust in most situations although they fail spectac- 
ularly on occasions. Reports on approximate Riemann solver failures and their respective 
corrections (usually a judicious addition of artificial dissipation) are abundant in the litera- 
ture |149|| . Motivated by the search for a robust and accurate approximate Riemann solver 



that avoids these common failures, Donat & Marquina E31 have extended to systems a 



numerical flux formula which was first proposed by Shu & Osher ||159| for scalar equations. 
In the scalar case and for characteristic wave speeds which do not change sign at the given 
numerical interface, Marquina's flux formula is identical to Roe's flux. Otherwise, the 



scheme switches to the more viscous, entropy satisfying local Lax-Friedrichs scheme [|159 
In the case of systems, the combination of Roe and local-Lax-Friedrichs solvers is carried 
out in each characteristic field after the local linearization and decoupling of the system of 



equations jK|. However, contrary to Roe's and other linearized methods, the extension of 



Marquina's method to systems is not based on any averaged intermediate state. 



Marti et al. have used this method in their simulations of relativistic jets ||107| , |108 |. 
The resulting numerical code has been successfully used to describe ultra-relativistic flows 
in both one and two spatial dimensions with great accuracy (a large set of test calculations 
using Marquina's Riemann solver can be found in Appendix II of [ |108|| ). Numerical exper- 
imentation in two dimensions confirms that the dissipation of the scheme is sufficient to 
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eliminate the carbuncle phenomenon ||149| , which appears in high Mach number relativistic 



jet simulations when using other standard solvers [42 1. Aloy et al. M have implemented 



Marquina's flux formula in their three dimensional relativistic hydrodynamic code GENE- 



SIS. Font et al. |59| have developed a 3D general relativistic hydro code where the matter 
equations are integrated in conservation form and fluxes are calculated with Marquina's 
formula. 

3.8 Symmetric TVD schemes with nonlinear numerical dissipa- 
tion 

The methods discussed in the previous subsections are all based on exact or approximate 
solutions of Riemann problems at cell interfaces in order to stabilize the discretization 
scheme across strong shocks. Another successful approach relies on the addition of non- 
linear dissipation terms to standard finite difference methods. The algorithm of Davis 



37J is based on such an approach. It can be interpreted as a Lax-Wendroff scheme with 
a conservative TVD dissipation term. The numerical dissipation term is local, free of 
problem dependent parameters and does not require any characteristic information. This 
last fact makes the algorithm extremely simple when applied to any hyperbolic system of 
conservation laws. 

A relativistic version of Davis' method has been used by Koide et al. |H], 80|, |126|| in 



2D and 3D simulations of relativistic magneto-hydrodynamic jets with moderate Lorentz 
factors. Although the results obtained are encouraging, the coarse grid zoning used in these 
simulations and the relative smallness of the beam flow Lorentz factor (4.56, beam speed 
~ 0.98c) does not allow for a comparison with Riemann-solver-based HRSC methods in 
the ultra-relativistic limit. 



4 OTHER DEVELOPMENTS 

4.1 Van Putten's approach 

Relying on a formulation of Maxwell's equations as a hyperbolic system in divergence form, 
van Putten [ 173 has devised a numerical method to solve the equations of relativistic ideal 



MHD in flat spacetime ||1 75| . Here we only discuss the basic principles of the method in 
one spatial dimension. In van Putten's approach, the state vector u and the fluxes F of 
the conservation laws are decomposed into a spatially constant mean (subscript 0) and a 
spatially dependent variational (subscript 1) part 

u{t,x) =u (t)+u 1 (t,x), F{t,x)=F (t)+F 1 (t,x). (43) 

The RMHD equations then become a system of evolution equations for the integrated 
variational parts ui*, which reads 

^ = 0, <«) 
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together with the conservation condition 

dt 

The quantities u x * are defined as 



. (45) 



/x 
Ui(t,y) dy. (46) 



They are continuous and standard methods can be used to integrate the system (|44|) . Van 
Putten uses a leapfrog method. 

The new state vector u(£, x) is then obtained from u-i*(t, x) by numerical differentiation. 
This process can lead to oscillations in the case of strong shocks and a smoothing algorithm 
should be applied. Details of this smoothing algorithm and of the numerical method in one 



and two spatial dimensions can be found in ||1 74j| together with results on a large variety 
of tests. 

Van Putten has applied his method to simulate relativistic hydrodynamic and magneto 
hydrodynamic jets with moderately flow Lorentz factors (< 4.25) ||1 76| , |179| . 



4.2 Relativistic SPH 

Besides finite volume schemes, another completely different method is widely used in as- 
trophysics for integrating the hydrodynamic equations. This method is Smoothed Particle 
Hydrodynamics, or SPH for short ( |]9~7| , [53| , |119| |). The fundamental idea of SPH is to 



represent a fluid by a Monte Carlo sampling of its mass elements. The motion and thermo- 
dynamics of these mass elements is then followed as they move under the influence of the 
hydrodynamics equations. Because of its Lagrangian nature there is no need within SPH 
for explicit integration of the continuity equation, but in some implementations of SPH 
this is done nevertheless for certain reasons. As both the equation of motion of the fluid 
and the energy equation involve continuous properties of the fluid and their derivatives, it 
is necessary to estimate these quantities from the positions, velocities and internal energies 
of the fluid elements, which can be thought of as particles moving with the flow. This 
is done by treating the particle positions as a finite set of interpolating points where the 
continuous fluid variables and their gradients are estimated by an appropriately weighted 
average over neighbouring particles. Hence, SPH is a free-Lagrange method, i.e., spatial 
gradients are evaluated without the use of a computational grid. 

A comprehensive discussion of SPH can be found in the reviews of Hernquist & Katz 
|76| , Benz |ll| and Monaghan [|118| , |119| . The non-relativistic SPH equations are briefly 



discussed in Section ( |9.5|) . The capabilities and limits of SPH are explored, e.g., in ||164| , 
167| , and the stability of the SPH algorithm is investigated in ||165| . 



The SPH equations for special relativistic flows have been first formulated by Monaghan 
118| . For such flows the SPH equations given in Section ( |9.5p can be taken over except 



that each SPH particle a carries v a baryons instead of mass m a ||118| , 29 . Hence, the rest 



mass of particle a is given by m a = mou a , where mo is the baryon rest mass (if the fluid is 
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made of baryons). Transforming the notation used in [^] to ours, the continuity equation, 
the momentum and the total energy equations for particle a are given by (unit of velocity 
is c) 

^r = - £ < w - - v ") • V ^ > ( 47 ) 



dt 

and 



r/S " ^ - l ''" + ^ + n ab ) ■ v a w ab , (« 



V- /Pa 



^a ^ /PaV a p b V b \ 



dt y V ^ ^ 



respectively Here, the summation is over all particles other than particle a, and d/dt 
denotes the Lagrangian time derivative. N = D/m is the baryon number density, 



the momentum per particle, and 



S = — = m hWv (50) 



T p 

— + m = m hW - — (51) 



the total energy per particle (all measured in the laboratory frame). The momentum 
density S = (S 1 , S 2 , S 3 ) T , the energy density r (measured in units of the rest mass energy 
density), and the specific enthalpy h are defined in Section fl2.1|) . II afe and fl a b are the SPH 
dissipation terms, and V W & denotes the gradient of the kernel W ab (see Section ( |9.5|) for 
more details). 



Special relativistic flow problems have been simulated with SPH by J89], 79, ^, |100| , |29 



|160| . Extensions of SPH capable of treating general relativistic flows have been considered 



by [79, |88], |160| . Concerning relativistic SPH codes the artificial viscosity is the most critical 



issue. It is required to handle shock waves properly, and ideally it should be predicted by 
a relativistic kinetic theory for the fluid. However, unlike its Newtonian analogue, the 
relativistic theory has not yet been developed to the degree required to achieve this. For 
Newtonian SPH Lattanzio et al. |)3| have shown that a viscosity quadratic in the velocity 
divergence is necessary in high Mach number flows. They proposed a form such that the 
viscous pressure could be simply added to the fluid pressure in the equation of motion 
and the energy equation. Because this simple form of the artificial viscosity has known 
limitations, they also proposed a more sophisticated form of the artificial viscosity terms, 
which leads to a modified equation of motion. This artificial viscosity works much better, 
but it cannot be generalized to the relativistic case in a consistent way. Utilizing an 



equation for the specific internal energy both Mann |99| and Laguna et al. |88| use such an 
inconsistent formulation. Their artificial visocity term is not included into the expression 
of the specific relativistic enthalpy. In a second approach, Mann |M|] allows for a time- 
dependent smoothing length and SPH particle mass, and further proposed a SPH variant 
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based on the total energy equation. Lahy |89[ and Siegler & Riffert |160|| use a consistent 
artificial viscosity pressure added to the fluid pressure. Siegler & Riffert 
formulated the hydrodynamic equations in conservation form. 



160 have also 



Monaghan ||120| incorporates concepts from Riemann solvers into SPH. For this reason 
he also proposes to use a total energy equation in SPH simulation instead of the commonly 
used internal energy equation, which would involve time derivatives of the Lorentz factor in 
the relativistic case. Chow & Monaghan [2]| have extended this concept and have proposed 
an SPH algorithm, which gives good results when simulating an ultra-relativistic gas. In 
both cases the intention was not to introduce Riemann solvers into the SPH algorithm, 
but to use them as a guide to improve the artificial viscosity required in SPH. 

as well as in its relativistic variant proposed by Eulerdink 
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In Roe's Riemann solver 
f50| , [5"Tf (see Section|3T4|), the numerical flux is computed by solving a locally linear system 
and depends on both the eigenvalues and (left and right) eigenvectors of the Jacobian 
matrix associated to the fluxes and on the jumps in the conserved physical variables (see 
Eqs. [26] and |27|). Monaghan ||120|| realized that an appropriate form of the dissipative terms 
Uab and Q a b for the interaction between particles a and b can be obtained by treating 
the particles as the equivalent of left and right states taken with reference to the line 
joining the particles. The quantity corresponding to the eigenvalues (wave propagation 
speeds) is an appropriate signal velocity v s i g (see below), and that equivalent to the jump 
across characteristics is a jump in the relevant physical variable. For the artificial viscosity 
tensor, II a b, Monaghan ||120|| assumes that the jump in velocity across characteristics can 
be replaced by the velocity difference between a and b along the line joining them. 

With these considerations in mind Chow & Monaghan [^] proposed for Ti^b i n the 
relativistic case the form 



n 



Kv sig (S* 



SJ)-j 



ab — T7 ) y oz ) 

^ ab 

when particles a and b are approaching, and II a & = otherwise. Here K = 0.5 is a 
dimensionless parameter, which is chosen to have the same value as in the non-relativistic 
case ||1 20[ . N a b = (N a + iV b )/2 is the average baryon number density, which has to be 
present in (|52"D , because the pressure terms in the summation of fl55|) have an extra density 



(53) 



in the denominator arising from the SPH interpolation. Furthermore, 

r ab 

J = I 1 

\ r ab\ 



is the unit vector from b to a, and 



where 



W* 



mt)hW*\ 



l-(v-j) 



(54) 



(55) 



Using instead of S (see Eq. [50]) the modified momentum S*, which involves the line of sight 
velocity v • j, guarantees that the viscous dissipation is positive definite [25|. 
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The dissipation term in the energy equation is derived in a similar way and is given by 

m 



IU = - *Wg-*»J , (56) 

^ ab 

if a and 6 are approaching, and Q a b = otherwise. Q a b involves the energy ?*, which is 
identical to f (see Eq. [[l|) except that W is replaced by W*. 

To determine the signal velocity Chow & Monaghan |^J (and Monaghan ||120|| in the 
non-relativistic case) start from the (local) eigenvalues, and hence the wave velocities (v ± 
c s )/(l ± vc s ) and v of one-dimensional relativistic hydrodynamic flows. Again considering 
particles a and b as the left and right states of a Riemann problem with respect to motions 
along the line joining the particles, the appropriate signal velocity is the speed of approach 
(as seen in the computing frame) of the signal sent from a towards b and that from b to a. 
This is the natural speed for the sharing of physical quantities, because when information 
about the two states meets it is time to construct a new state. This speed of approach 
should be used when determining the size of the time step by the Courant condition (for 
further details see [p9J1 ). 



Chow & Monaghan |^9[ have demonstrated the performance of their Riemann problem 



guided relativistic SPH algorithm by calculating several shock tube problems involving 
ultra-relativistic speeds up to v — 0.9999. The algorithm gives good results, but finite 
volume schemes based on Riemann solvers give more accurate results and can handle even 
larger speeds (see Section^). 

4.3 Relativistic beam scheme 



Sanders & Prendergast |155|| proposed an explicit scheme to solve the equilibrium limit 



of the non-relativistic Boltzmann equation, i.e., the Euler equations of Newtonian fluid 
dynamics. In their so-called beam scheme the Maxwellian velocity distribution function 
is approximated by several Dirac delta functions or discrete beams of particles in each 
computational cell, which reproduce the appropriate moments of the distribution function. 
The beams transport mass, momentum and energy into adjacent cells, and their motion 
is followed to first-order accuracy. The new (i.e., time advanced) macroscopic moments 
of the distribution function are used to determine the new local non-relativistic Maxwell 
distribution in each cell. The entire process is then repeated for the next time step. The 
CFL stability condition requires that no beam of gas travels farther than one cell in one 
time step. This beam scheme, although being a particle method derived from a microscopic 
kinetic description, has all the desirable properties of modern characteristic-based wave 
propagating methods based on a macroscopic continuum description. 

The non-relativistic scheme of Sanders & Prendergast [|155|| has been extended to rel- 
ativistic flows by Yang et ai. |191[ | . They replaced the Maxwellian distribution function 



by its relativistic analogue, i.e., by the more complex Jiittner distribution function, which 
involves modified Bessel functions. For three-dimensional flows the Jiittner distribution 
function is approximated by seven delta functions or discrete beams of particles, which 
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can viewed as dividing the particles in each cell into seven distinct groups. In the local 
rest frame of the cell these seven groups represent particles at rest and particles moving in 
±x, ±y and ±z directions, respectively. 

Yang et al. ||191|1 show that the integration scheme for the beams can be cast in the 
form of an upwind conservation scheme in terms of numerical fluxes. They further show 
that the beam scheme not only splits the state vector but also the flux vectors, and has 
some entropy-satisfying mechanism embedded as compared with approximate relativistic 
Riemann solver [11, |157|| based on Roe's method ||151 |. The simplest relativistic beam 



scheme is only first-order accurate in space, but can be extended to higher-order accu- 
racy in a straightforward manner. Yang et al. consider three high-order accurate variants 
(TVD2, EN02, EN03) generalizing their approach developed in ||189| , |190|| for Newtonian 



gas dynamics, which is based on the essentially non-oscillatory (ENO) piecewise polynomial 
reconstruction scheme of Harten et al. [[73] . 



Yang et al. [ p.91| ] present several numerical experiments including relativistic one-dimen- 
sional shock tube flows and the simulation of relativistic two-dimensional Kelvin-Helmholtz 
instabilities. The shock tube experiments consist of a mildly relativistic shock tube, rel- 
ativistic shock heating of a cold flow, the relativistic blast wave interaction of Woodward 
& Colella [|186|| (see Section |6.2.3j ), and the perturbed relativistic shock tube flow of Shu & 



Osher 159 



5 SUMMARY OF METHODS 

This Section contains a summary of all the methods reviewed in the two preceding sections 
as well as several FCT and artificial viscosity codes. The main characteristic of the codes 
(dissipation algorithm, spatial and temporal orders of accuracy, reconstruction techniques) 
are listed in two tables (Table |2| for HRSC codes; Table ||] for other approaches). 

6 TEST BENCH 

6.1 Relativistic shock heating in planar, cylindrical and spherical 
geometry 

Shock heating of a cold fluid in planar, cylindrical or spherical geometry has been used 
since the early developments of numerical relativistic hydrodynamics as a test case for 
hydrodynamic codes, because it has an analytical solution ([[L7[] in planar symmetry; ||108 



in cylindrical and spherical symmetry), and because it involves the propagation of a strong 
relativistic shock wave. 

In planar geometry, an initially homogeneous, cold (i.e.,e ~ 0) gas with coordinate 
velocity v i and Lorentz factor W\ is supposed to hit a wall, while in the case of cylindrical 
and spherical geometry the gas flow converges towards the axis or the center of symmetry. 
In all three cases the reflection causes compression and heating of the gas as kinetic energy 
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Table 2: High-resolution shock-capturing methods. All the codes rely on a conservation 



form of the RHD equations with the exception of ref. ||182 



Code 



Basic characteristics 



Roe type-1 [104, 153, F>9| 



Roe-Eulderink 



HLL-1 [157| 



LCA-phm [103 1 



LCA-eno @ 



rPPM [106| 



Falle-Komissarov [55 
MFF-ppm JHJ|, | 



MFF-eno/phm |§ 



MFF-1 



Flux split |>G 
sTVD §1 



rGlimm |182| 
rBS Il9l| 



Riemann solver of Roe type with arithmetic averaging; monotonicity 
preserving, linear reconstruction of primitive variables; 2nd order time 
stepping ([104, 153 1: predictor-corrector; |5J|: standard scheme) 

Linearized Riemann solver based on Roe averaging; 2nd order 
accuracy in space and timere 

Harten-Lax-van Leer approximate Riemann solver; monotonic linear 
reconstruction of conserved/primitive variables; 2nd order accuracy 
in space and time 

Local linearization and decoupling of the system; PHM reconstruction 
of characteristic fluxes; 3rd order TVD preserving RK method for 
time stepping 

Local linearization and decoupling of the system; high order ENO 
reconstruction of characteristic split fluxes; high order TVD 
preserving RK methods for time stepping 

Exact (ideal gas) Riemann solver; PPM reconstruction of primitive 
variables; 2nd order accuracy in time by averaging states in the 
domain of dependence of zone interfaces 

Approximate Riemann solver based on local linearizations of the 
RHD equations in primitive form; monotonic linear reconstruction 
of p, p and u l ] 2nd order predictor-corrector time stepping 

Marquina flux formula for numerical flux computation; PPM 
reconstruction of primitive variables; 2nd and 3rd order TVD 
preserving RK methods for time stepping 

Marquina flux formula for numerical flux computation; upwind biased 
ENO/PHM reconstruction of characteristic fluxes; 2nd and 3rd order 
TVD preserving RK methods for time stepping 

Marquina flux formula for numerical flux computation; monotonic 
linear reconstruction of primitive variables; standard 2nd order 
finite difference algorithms for time stepping 

TVD flux-split second order method 

Davis (1984) symmetric TVD scheme with nonlinear numerical 

dissipation; 2nd order accuracy in space and time 

Glimm's method applied to RHD equations in primitive form; 1st 
order accuracy in space and time 

Relativistic beam scheme solving equilibrium limit of relativistic 
Boltzmann equation; distribution function approximated by discrete 
beams of particles reproducing appropriate moments; 1st and 2nd 
order TVD, 2nd and 3rd order ENO schemes 
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Table 3: Code characteristics 



Code 



Basic characteristics 



Artificial viscosity 



AV-mono [gfl |75J |110| 



cAV-implicit 02 



Non-conservative formulation of the RHD equations (transport 
differencing, internal energy equation); artificial viscosity 
extra term in the momentum flux; monotonic 2nd order 
transport differencing; explicit time stepping 

Non-conservative formulation of the RHD equations; internal 
energy equation; consistent formulation of artificial viscosity; 
adaptive mesh and implicit time stepping 



Flux corrected transport 



FCT-lw 44 



SHASTA-c |157], |38|, |3 



Non-conservative formulation of the RHD equations (transport 
differencing, equation for phW); explicit 2nd order 
Lax-Wendroff scheme with FCT algorithm 
FCT algorithm based on SHASTA f9); advection of 
conserved variables 



van Putten's approach 



van Putten fll7 



Ideal RMHD equations in constraint-free, divergence form; 
evolution of integrated variational parts of conserved quantities; 
smoothing algorithm in numerical differentiation step; 
leap-frog method for time stepping 



Smooth particle hydrodynamics 



SPH-AV-0 @(SPH0), 



SPH-AV-1 H(SPH1) 



SPH-AV-c ||(SPH2) 



SPH-cAV-c [160 1 



SPH-RS-c 2£ 



Specific internal energy equation; artificial viscosity extra terms 
in momentum and energy equations; 2nd order time stepping 
(|99[|: predictor-corrector; [p8[ : RK method) 

Time derivatives in SPH equations include variations in smoothing 

length and mass per particle; Lorentz factor terms treated more 

consistently; otherwise same as SPH-AV-0 

Total energy equation; otherwise same as SPH-AV-1 

RHD equations in conservation form; consistent formulation of 

artificial viscosity 

RHD equations in conservation form; dissipation terms constructed 
in analogy to terms in Ricmann solver based methods 
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is converted into internal energy. This occurs in a shock wave, which propagates upstream. 
Behind the shock the gas is at rest (v 2 = 0). Due to conservation of energy across the 
shock the gas has a specific internal energy given by 

e 2 = W 1 -l. (57) 

The compression ratio of shocked and unshocked gas, a, follows from 

a = ^±± + ^-e 2 , (58) 

7 — 1 7 — 1 

where 7 is the adiabatic index of the equation of state. The shock velocity is given by 

(7-l)Wi|vi| 



U 



Wi + 1 



(59) 



In the unshocked region (r e [V s t, 00 [) the pressureless gas flow is self-similar and has a 
density distribution given by 



p(t,r) 




(60) 



where a = 0, 1, 2 for planar, cylindrical or spherical geometry, and where po is the density 
of the inflowing gas at infinity (see Fig. ^) . 

In the Newtoninan case the compression ratio a of shocked and unshocked gas cannot 
exceed a value of cr max = (7 + l)/(7 — 1) independently of the inflow velocity. This is 
different for relativistic flows, where a grows linearly with the flow Lorentz factor and 
becomes infinite as the inflowing gas velocity approaches to speed of light. 

The maximum flow Lorentz factor achievable for a hydrodynamic code with acceptable 
errors in the compression ratio a is a measure of the code's quality. Table|] contains a 
summary of the results obtained for the shock heating test by various authors. 

Explicit finite-difference techniques based on a non-conservative formulation of the 
hydrodynamic equations and on non-consistent artificial viscosity |27, |75| are able to handle 
flow Lorentz factors up to ~ 10 with moderately large errors (o" err0 r ~ 1 — 3%) at best 
185], |110|| . Norman & Winkler |128|| got very good results (cr crror ~ 0.01% for a flow 
Lorentz factor of 10 using consistent artificial viscosity terms and an implicit adaptive- 
mesh method. 

The performance of explicit codes improved significantly when numerical methods based 
on Riemann solvers were introduced [|104| , |103| , |50| , |157| , |5T| , |106| , [55| . For some of these codes 
the maximum flow Lorentz factors is only limited by the precision by which numbers are 
represented on the computer used for the simulation |y], |182| , |2j . 



Schneider et al. |157|| have compared the accuracy of a code based on the relativistic 
HLL Riemann solver with different versions of relativistic FCT codes for inflow Lorentz 
factors in the range 1.6 to 50. They found that the error in a was reduced by a factor of 
two when using HLL. 
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Relativistic shock heating 
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Figure 3: Schematic solution of the shock heating problem in spherical geometry. The 
initial state consists of a spherically symmetric flow of cold (p = 0) gas of unit rest mass 
density having a coordinate velocity v\ = —1 everywhere. A shock is generated at the 
center of the sphere, which propagates upstream with constant speed. The post-shock 
state is constant and at rest. The pre-shock state, where the flow is self-similar, has a 
density which varies as p = (1 + t/r) 2 with time t and radius r. 
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Within SPH methods, Chow & Monaghan [29| have obtained results comparable to 
those of HRSC methods (cr error < 2 10~ 3 ) for flow Lorentz factors up to 70, using a rela- 



tivistic SPH code with Riemann solver guided dissipation. Sieglert & Riffert [160] have 
succeded in reproducing the post-shock state accurately for inflow Lorentz factors of 1000 
with a code based on a consistent formulation of artificial viscosity. However, the disspa- 
tion introduced by SPH methods at the shock transition is very large (10-12 particles in 
the code of ref. ||160[| ; 20-24 in the code of ref. p9| ) compared with the typical dissipation 



of HRSC methods (see below). 

The performance of a HRSC method based on a relativistic Riemann solver is illustrated 
by means of a MPEG movie (final frame of movie is displayed in Fig. 0) for the planar 
shock heating problem for an inflow velocity V\ = —0.99999c (W\ ~ 223). These results 
are obtained with the relativistic PPM code of |106|| , which uses an exact Riemann solver 



based on the procedure described in Section 2.3. 



The shock wave is resolved by three zones and there are no post-shock numerical 
oscillations. The density increases by a factor p» 900 across the shock. Near x = 
the density distribution slightly undershoots the analytical solution (by « 8%) due to 
the numerical effect of wall heating. The profiles obtained for other inflow velocities are 
qualitatively similar. The mean relative error of the compression ratio 0" errO r < 10~ 3 , and, 
in agreement with other codes based on a Riemann solver, the accuracy of the results does 
not exhibit any significant dependence on the Lorentz factor of the inflowing gas. 

Some authors have considered the problem of shock heating in cylindrical or spherical 
geometry using adapted coordinates to test the numerical treatment of geometrical factors 
|153| , |108| , |182j ]. Aloy et al. [[| have considered the spherically symmetric shock heating 



problem in 3D Cartesian coordinates as a test case for both the directional splitting and 
the symmetry properties of their code GENESIS. The code is able to handle this test up 
to inflow Lorentz factors of the order of 700. 

In the shock reflection test conventional schemes often give numerical approximations 
which exhibit a consistent 0(1) error for the density and internal energy in a few cells near 
the reflecting wall. This 'overheating', as it is known in classical hydrodynamics ||127| , is a 
numerical artifact which is considerably reduced when Marquina's scheme is used [ O . 

6.2 Propagation of relativistic blast waves 

Riemann problems with large initial pressure jumps produce blast waves with dense shells 
of material propagating at relativistic speeds (see Fig. [5]). For appropriate initial conditions, 
both the speed of the leading shock front and the velocity of the shell material approach 
the speed of light producing very narrow structures. The accurate description of these 
thin, relativistic shells involving large density contrasts is a challenge for any numerical 
code. Some particular blast wave problems have become standard numerical tests. Here we 
consider the two most common of these tests. The initial conditions are given in Table^. 
Problem 1 was a demanding problem for relativistic hydrodynamic codes in the mid 



eighties |27|, [75|], while Problem 2 is a challenge even for today's state-of-the-art codes. The 



analytical solution of both problems can be obtained with program RIEMANN (see Sect. 
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Table 4: Summary of relativistic shock heating test calculations by various authors in planar (a = 
0), cylindrical (a = 1), and spherical (a = 2) geometry. VF max and <7 error are the maximum inflow 
Lorentz factor and compression ratio error extracted from tables and figures of the corresponding 
reference. W max should only be considered as indicative of the maximum Lorentz factor achievable 
by every method. Methods are described in Sects. || and ^| and their basic properties summarized 
in Sect. | (Tables ID). 



References 



a Method 



W n 



<Jc 



Centrella & Wilson (1984) [|7]] 
Hawley et al. (1984) J7| 
Norman & Winkler (1986) pg] 
McAbee et al (1989) pTJfl 
Marti et al (1991) M 
Marquina et al (1992) flTUB" 
Eulderink (1993) [0 
Schneider et al. (1993) PI 



Dolezal & Wong (1995) g| 
Marti & Miiller (1996) |T06f 
Falle & Komissarov (1996) 
Romero et al. (1996) 



Marti et al. (1997) [ |10l 
Chow & Monaghan (1997) 
Wen et al. (1997) p 
Donat et al. (1998) 
Aloy et al. (1999) [ 



Sieglert & Riffert (1999) |60| 
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AV-mono 


4.12 


w 10 





cAV-implicit 


10.0 


0.01 





AV-mono 


10.0 


2.6 





Roe type-1 


23 


0.2 





LCA-phm 


70 


0.1 





Roe-Eulderink 


625 


< 0.1 a 





HLL-1 


10 6 


0.2 b 





SHASTA-c 


10 6 


0.5 b 





LCA-eno 


7.0 10 5 


< 0.1 a 





rPPM 


224 


0.03 





Falle- Komissarov 


224 


< 0.1 a 


2 


Roe type-1 


2236 


2.2 


1 


MFF-ppm 


70 


1.0 





SPH-RS-c 


70 


0.2 


2 


rGlimm 


224 


1Q -9 





MFF-eno 


224 


< 0.1 a 





MFF-ppm 


2.4 10 5 


3.5 C 





SPH-cAV-c 


1000 


< 0.1 a 



a Estimated from figures. 

b For !T max = 50. 

c Including points at shock transition. 
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Figure 4: MPEG movie (only final frame of movie is displayed f) showing the evolution of the 
density distribution for the shock heating problem with an inflow velocity v\ = —0.99999c in 
Cartesian coordinates. The reflecting wall is located at x = 0. The adiabatic index of the gas is 
4/3. For numerical reasons, the specific internal energy of the inflowing cold gas is set to a small 
finite value (ei = 10 -7 W\). The final frame of the movie also shows the analytical solution (blue 
lines). The simulation has been performed on an equidistant grid of 100 zones. 
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Generation and propagation of 
relativistic blast waves 
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Figure 5: Generation and propagation of a relativistic blast wave (schematic). The large pressure 
jump at a discontinuity initially located at r = 05 gives rise to a blast wave and a dense shell of 
material propagating at relativistic speeds. For appropriate initial conditions both the speed of 
the leading shock front and the velocity of the shell approach the speed of light producing very 
narrow structures. 
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Table 5: Initial data (pressure p, density p, velocity v) for two common relativistic blast wave test 
problems. The decay of the initial discontinuity leads to a shock wave (velocity v s hock) compression 
ratio cr s hock) and the formation of a dense shell (velocity i> s heii> time-dependent width u> s heii) both 
propagating to the right. The gas is assumed to be ideal with an adiabatic index 7 = 5/3. 



Problem 1 Problem 2 

Left Right Left Right 



V 


13.33 0.00 


1000.00 0.01 


p 


10.00 1.00 


1.00 1.00 


V 


0.00 0.00 


0.00 0.00 


^shell 


0.72 


0.960 


Wshell 


O.llt 


0.026 1 


^shock 


0.83 


0.986 


""shock 


5.07 


10.75 



6.2.1 Problem 1 

In Problem 1 , the decay of the initial discontinuity gives rise to a dense shell of matter 
with velocity w s h e ii = 0.72 (W B heii = 1-38) propagating to the right. The shell trailing a 
shock wave of speed f s hock = 0.83 increases its width, u> s heii, according to u> s heii = O.llt, 
i.e., at time t = 0.4 the shell covers about 4% of the grid (0 < x < 1). Tables^?] give 
a summary of the references where this test was considered for non-HRSC and HRSC 
methods, respectively. 

Using artificial viscosity techniques Centrella & Wilson |27| were able to reproduce the 



analytical solution with a 7% overshoot in f s h e ii, whereas Hawley et al. |75[ got a 16% error 
in the shell density. 

The results obtained with early relativistic SPH codes |M| were affected by systematic 
errors in the rarefaction wave and the constant states, large amplitud spikes at the contact 
discontinuity and large smearing. Smaller systematic errors and spikes are obtained with 



Laguna et al. 's (1993) code [ 88| . This code also leads to a large overshoot in the shell's 
density. Much cleaner states are obtained with the methods of Chow & Monaghan (1997) 
and Siegler & Riffert (1999) [|160|| , both based on conservative formulations of the 



SPH equations. In the case of Chow & Monaghan's (1997) method [p29 , the spikes at the 



contact discontinuity dissapear but at the cost of an excessive smearing. Shock profiles with 



32 



Table 6: Summary of references where the blast wave Problem 1 (defined in Tabled) has been 
considered in ID, 2D and 3D, respectively. Methods are described in Sects. |3| and || and their 
basic properties summarized in Sect. || (Tables |2]||). 



References 



Dim. Method 



Comments 



Centrella & Wilson (1984) | 
Hawley et al. (1984) || 
Dubai (1991) a || 



Mann (1991) || 

Laguna et al. (1993) |8 
van Putten (1993) 6 |7 



Schneider et al. (1993) p7j 



Chow & Monaghan (1997) || 



Siegler & Riffert (1999) JT6fJ | 



ID AV-mono 

ID AV-mono 

ID FCT-lw 

ID SPH-AV-0,1,2 



ID 


SPH-AV-0 


ID 


van Putten 


ID 


SHASTA-c 


ID 


SPH-RS-c 


ID 


SPH-cAV-c 



Stable profiles without oscillations 

Velocity overestimated by 7% 

Stable profiles without oscillations 

Psheii overestimated by 16% 

10-12 zones at the CD 

Velocity overestimated by 4.5% 

Sistematic errors in the rarefaction 

wave and the constant states 

Large amplitud spikes at the CD 

Excessive smearing at the shell 

Large amplitud spikes at the CD 

Psheii overestimated by 5% 

Stable profiles 

Excessive smearing, specially at the 

CD (« 50 zones) 

Non monotonic intermediate states 

p s h e ii underestimated by 10% with 

200 zones 

Stable profiles without spikes 

Excessive smearing at the CD and 

at the shock 

Correct constant states 

Large amplitud spikes at the CD 

Excessive smearing at the shock 

transition (^ 20 zones) 



a For a Riemann problem with slightly different initial conditions. 

b For a Riemann problem with slightly different initial conditions including a nonzero 



transverse magnetic field. 
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Table 7: Summary of references where the blast wave Problem 1 (defined in Tableg) has been 
considered in ID, 2D and 3D, respectively. Methods are described in Sects. |3] and [| and their 
basic properties summarized in Sect. || (Tables |2|,y). 



References 



Dim. Method 



Comments 



Eulderink (1993) |(J 
Schneider et al (1993) |5 
Marti & Miiller (1996) flIP5 | 
Marti et al (1997) |D| 



Wen et al. (1997) [ |T82j 
Yang et al (1997) |9|] 
Donat et al. (1998) [[TJ 

Aloy et al. (1999) § 
Font et al. (1999) pf 



ID 


Roe-Eulderi 


ID 


HLL-1 


ID 


rPPM 


ID, 2D 


MFF-ppm 


ID 


rGlimm 


ID 


rBS 


ID 


MFF-eno 


3D 


MFF-ppm 


ID, 3D 


MFF-1 


ID, 3D 


Roe type-1 


ID, 3D 


Flux split 



Correct p s heii with 500 zones 

4 zones at the CD 

Pshdi underestimated by 10% 

with 200 zones 

Correct p s heii with 400 zones 

6 zones at the CD 

Correct p s heii with 400 zones 

6 zones at the CD 

No difussion at discontinuities 

Stable profiles 

Correct p s heii with 400 zones 

8 zones at the CD 

Correct p s heii with 400 zones 
12-14 zones at the CD 
Correct p s heii with 400 zones 
12-14 zones at the CD 
Psheii overestimated by 5% 
8 zones at the CD 



a All the methods produce stable profiles without numerical oscillations, 
correspond to ID cases. 



Comments 
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relativistic SPH codes are more smeared out than with HRSC methods covering typically 
more than 10 zones. 

Van Putten has considered a similar initial value problem with somewhat more extreme 
conditions (f shell ~ 0.82c, cr s h oc k ~ 5.1) and with a transversal magnetic field. For suitable 
choices of the smoothing parameters his results are accurate and stable, although discon- 
tinuities appear to be more smeared than with typical HRSC methods (6-7 zones for the 
strong shock wave; ~ 50 zones for the contact discontinuity). 

MPEG movie 2 (final frame of movie is displayed in Fig. |6|) shows the Problem 1 blast 
wave evolution obtained with a modern HRSC method (the relativistic PPM method in- 



troduced in Section pTT]) . The grid has 400 equidistant zones, and the relativistic shell is 
resolved by 16 zones. Because of both the high order accuracy of the method in smooth 
regions and its small numerical diffusion (the shock is resolved with 4-5 zones only) the 
density of the shell is accurately computed (errors less than 0.1%). Other codes based 
on relativistic Riemann solvers |H| give similar results (see Table |7j). The relativistic 



HLL method |157|| underestimates the density in the shell by about 10% in a 200 zone 
calculation. 

6.2.2 Problem 2 



Problem 2 was first considered by Norman & Winkler ||1 28|| . The flow pattern is similar 
to that of Problem 1, but more extreme. Relativistic effects reduce the post -shock state 
to a thin dense shell with a width of only about 1% of the grid length at t = 0.4. The 
fluid in the shell moves with t> s h e ii = 0.960 (i.e., W^heii = 3.6), while the leading shock front 
propagates with a velocity t> s hock = 0.986 (i.e., W^hock = 6.0). The jump in density in the 



shell reaches a value of 10.6. Norman & Winkler ||128| obtained very good results with 



an adaptive grid of 400 zones using an implicit hydro-code with artificial viscosity. Their 
adaptive grid algorithm placed 140 zones of the available 400 zones within the blast wave 
thereby accurately capturing all features of the solution. 

Several HRSC methods based on relativistic Riemann solvers have used Problem 2 as 
a standard test |0|, [10|, [T06], g5] $M M (see Table |). 



MPEG movie 3 (final frame of movie is displayed in Fig. [?J) shows the Problem 2 blast 



wave evolution obtained with the relativistic PPM method introduced in Section 34.) on 
a grid of 2000 equidistant zones. At this resolution the relativistic PPM code obtains a 
converged solution. The method of Falle & Komissarov J5l| requires a seven-level adaptive 
grid calculation to achieve the same the finest grid spacing corresponding to a grid of 3200 
zones. As their code is free of numerical diffusion and dispersion, Wen et ai. ||182| are able 
to handle this problem with high accuracy (see Fig[|). At lower resolution (400 zones) the 
relativistic PPM method only reaches 69% of the theoretical shock compression value (54% 



in case of the second-order accurate upwind method of Falle & Komissarov |55| ; 60% with 



the code of Donat et ai. [|42|] ) . 

Chow & Monaghan |29| have considered Problem 2 to test their relativistic SPH code. 
Besides a 15% overshoot in the shell's density, the code produces a non-causal blast wave 
propagation speed (i.e.,t> s h oc k > !)• 
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Figure 6: MPEG movie {only final frame of movie is displayed!) showing the evolution of the 
density distribution for the relativistic blast wave Problem 1 defined in Table |5[ The final frame 
of the movie also shows the analytical solution (blue lines). The simulation has been performed 
with relativistic PPM on an equidistant grid of 400 zones. 



36 



Table 8: Summary of references where the blast wave Problem 2 (defined in Table ||) has been 
considered. Methods are described in Sects. |3| and || and their basic properties summarized in 
Sect. | (Tables|2]||). 



References 


Method 


°7°cxact 


Norman & Winkler (1986) |128 


cAV-implicit 


1.00 


Dubai (1991) |4| a 


FCT-lw 


0.80 


Marti et al. (1991) ||104| 


Roe type-1 


0.53 


Marquina et al. (1992) pTJ3|] 


LCA-phm 


0.64 


Marti &Muller (1996) ||106| 


rPPM 


0.68 


Falle & Komissarov (1996) |5§ 


Falle- Komissarov 


0.47 


Wenet al. (1997) ||182|| 


rGlimm 


1.00 


Chow & Monaghan (1997) || 


SPH-RS-c 


1.16 6 


Donat et al. (1998) [0 


MFF-phm 


0.60 



a For a Riemann problem with slightly different initial conditions. 
b At t = 0.15. 
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Figure 7: MPEG movie (only final frame of movie is displayed !) showing the evolution of the 
density distribution for the relativistic blast wave Problem 2 defined in Table |5[ The final frame 
of the movie also shows the analytical solution (blue lines). The simulation has been performed 
with relativistic PPM on an equidistant grid of 2000 zones. 
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Relativistic shock tube 1 



Relativistic shock tube 2 
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Figure 8: Results from [182] for the relativistic blast wave Problems 1 (left column) and 2 (right 
column), respectively. Relativistic Glimm's method is only used in regions with steep gradients. 
Standard finite difference schemes are applied in the smooth remaining part of the computational 
domain. In the above plots, Lax and LW stand respectively for Lax and Lax-Wendroff methods; 
G refers to pure Glimm's method. 
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6.2.3 Collision of two relativistic blast waves 

The collision of two strong blast waves was used by Woodward & Colella [|186|| to compare 
the performance of several numerical methods in classical hydrodynamics. In the rela- 



tivistic case, Yang et al. ||191|| considered this problem to test the high-order extensions of 



the relativistic beam scheme, whereas Marti & Miiller ||106| used it to evaluate the perfor- 
mance of their relativistic PPM code. In this last case, the original boundary conditions 
were changed (from reflecting to outflow) to avoid the reflection and subsequent interaction 
of rarefaction waves allowing for a comparison with an analytical solution. In the following 
we summarize the results on this test obtained by Marti & Miiller in [|106|| . 

The initial data corresponding to this test, consisting in three constant states with large 
pressure jumps at the discontinuities separating the states (at x = 0.1 and x = 0.9), as 
well as the properties of the blast waves created by the decay of the initial discontinuities, 
are listed in Table ^. The propagation velocity of the two blast waves is slower than in 
the Newtonian case, but very close to the speed of light (0.9776 and —0.9274 for the shock 
wave propagating to the right and left, respectively). Hence, the shock interaction occurs 
later (at t = 0.420) than in the Newtonian problem (at t = 0.028). Top panel in Fig. ^ 
shows four snapshots of the density distribution including the moment of the collision of 
the blast waves at t = 0.420 and x = 0.5106. At the time of collision the two shells have a 
width of Ax = 0.008 (left shell) and Ax = 0.019 (right shell), respectively, i.e., the whole 
interaction takes place in a very thin region (about 10 times smaller than in the Newtonian 
case where Ax ~ 0.2). 

The collision gives rise to a narrow region of very high density (see lower panel of Fig. ^) 
bounded by two shocks moving at speeds 0.088 (shock at the left) and 0.703 (shock at the 
right) and large compression ratios (7.26 and 12.06, respectively) well above the classical 
limit for strong shocks (6.0 for 7 = 1.4). The solution just described applies until t = 0.430 
when the next interaction takes place. 

The complete analytical solution before and after the collision up to time t = 0.430 can 
be obtained following Appendix II in | |106[ . 



MPEG movie 4 (final frame of movie is displayed in Fig. \ffi) shows the evolution of the 
density up to the time of shock collision at t — 0.4200. The movie was obtained with the 



relativistic PPM code of Marti & Miiller [IOC]. The presence of very narrow structures 



involving large density jumps requires very fine zoning to resolve the states properly. For 
the movie a grid of 4000 equidistant zones was used. The relative error in the density 
of the left (right) shell is always less than 2.0% (0.6%), and is about 1.0% (0.5%) at the 
moment of shock collision. Profiles obtained with the relativistic Godunov method (first- 
order accurate, not shown) show relative errors in the density of the left (right) shell of 
about 50% (16%) at t — 0.20. The errors drop only slightly to about 40% (5%) at the time 
of collision (t = 0.420). 

MPEG movie 5 (final frame of movie is displayed in Fig. [Tip shows the numerical solu- 
tion after the interaction has occurred. Compared to MPEG movie 4 (final frame of movie 
is displayed in Fig. |I7J) a very different scaling for the x-axis had to be used to display 
the narrow dense new states produced by the interaction. Obviously, the relativistic PPM 
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Interaction of Relativistic Blast Waves 
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Figure 9: The top panel shows a sequence of snapshots of the density profile for the colliding 
relativistic blast wave problem up to the moment when the waves begin to interact. The density 
profile of the new states produced by the interaction of the two waves is shown in the bottom 
panel (note the change in scale on both axes with respect to the top panel). 
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Figure 10: MPEG movie (only final frame of movie is displayed!) showing the evolution of 
the density distribution for the colliding relativistic blast wave problem up to the interaction of 
the waves. The final frame of the movie also shows the analytical solution (blue lines). The 
computation has been performed with relativistic PPM on an equidistant grid of 4000 zones. 
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Figure 11: MPEG movie (only final frame of movie is displayed!) showing the evolution of the 
density distribution for the colliding relativistic blast wave problem around the time of inter- 
action of the waves at an enlarged spatial scale. The final frame of the movie also shows the 
analytical solution (blue lines). The computation has been performed with relativistic PPM on 
an equidistant grid of 4000 zones. 
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Table 9: Initial data (pressure p, density p, velocity v) for the two relativistic blast wave collision 
test problem. The decay of the initial discontinuities (at x = 0.1 and x = 0.9) leads to the 
formation of two shock waves (velocities w s hock, compression ratios <7 s hock) leading two dense 
shells (velocities t> s h e ii> time-dependent widths u? s heii) moving in opposite directions. The gas is 
assumed to be ideal with an adiabatic index 7 = 1.4. 



Left Middle Right 



p 


1000.00 


0.01 100.00 


p 


1.00 


1.00 1.00 


V 


0.00 


0.00 0.00 


^shell 


0.957 


-0.882 


Wshell 


0.021 1 


0.045 1 


^shock 


0.978 


-0.927 


C'shock 


14.39 


9.72 



code resolves the structure of the collision region satisfactorily well, the maximum relative 
error in the density distribution being less than 2.0%. When using the first-order accurate 
Godunov method instead, the new states are strongly smeared out and the positions of the 
leading shocks are wrong. 

7 APPLICATIONS 

7.1 Astrophysical jets 

The most compelling case for a special relativistic phenomenon are the ubiquitous jets 
in extragalactic radio sources associated with active galactic nuclei. In the commonly 
accepted standard model |TD[, flow velocities as large as 99% (in some cases even beyond) 



of the speed of light are required to explain the apparent superluminal motion observed 
in many of these sources. Models which have been proposed to explain the formation of 
relativistic jets, involve accretion onto a compact central object, such as a neutron star or 



stellar mass black hole in the galactic micro-quasars GRS 1915+105 Hi 16(1 an d GRO J1655 



40 1 169 1 , or a rotating super massive black hole in an active galactic nucleus, which is fed 
by interstellar gas and gas from tidally disrupted stars. 

Inferred jet velocities close to the speed of light suggest that jets are formed within 
a few gravitational radii of the event horizon of the black hole. Moreover, very-long- 

44 



baseline interferometric (VLBI) radio observations reveal that jets are already collimated 
at subparsec scales. Current theoretical models assume that accretion disks are the source 
of the bipolar outflows which are further collimated and accelerated via MHD processes 
(see, e.g., [15|]). There is a large number of parameters which are potentially important for 



jet powering: the black hole mass and spin, the accretion rate and the type of accretion 
disk, the properties of the magnetic field and of the environment. 

At parsec scales the jets, observed via their synchrotron and inverse Compton emission 
at radio frequencies with VLBI imaging, appear to be highly collimated with a bright 
spot (the core) at one end of the jet and a series of components which separate from 



the core, sometimes at superluminal speeds. In the standard model |16| , these speeds are 
interpreted as a consequence of relativistic bulk motions in jets propagating at small angles 
to the line of sight with Lorentz factors up to 20 or more. Moving components in these 
jets, usually preceded by outbursts in emission at radio wavelengths, are interpreted in 
terms of traveling shock waves. 

Finally, the morphology and dynamics of jets at kiloparsec scales are dominated by 
the interaction of the jet with the surrounding extragalactic medium, the jet power being 
responsible for dichotomic morphologies (the so called Fanaroff-Riley I and II classes p6| , 



FRI and FRII, respectively). Whereas current models [T3|, |K| interpret FRI morphologies 
as the result of a smooth deceleration from relativistic to non-relativistic, transonic speeds 
on kpc scales due to a slower shear layer, flux asymmetries between jets and counter-jets 
in the most powerful radio galaxies (FRII) and quasars indicate that relativistic motion 
extends up to kpc scales in these sources, although with smaller values of the overall bulk 
speeds |f2U . 



Although MHD and general relativistic effects seem to be crucial for a successful launch 
of the jet (for a review see, e.g., [E2|), purely hydrodynamic, special relativistic simula- 
tions are adequate to study the morphology and dynamics of relativistic jets at distances 
sufficiently far from the central compact object (i.e., at parsec scales and beyond). The 
development of relativistic hydrodynamic codes based on HRSC techniques (see Sections^ 
and H) has triggered the numerical simulation of relativistic jets at parsec and kiloparsec 
scales. 

At kiloparsec scales, the implications of relativistic flow speeds and/or relativistic in- 
ternal energies for the morphology and dynamics of jets have been the subject of a number 
of papers in recent years ||109| , [46] , |107[ |108| , J35f . Beams with large internal energies show 



little internal structure and relatively smooth cocoons allowing the terminal shock (the 
hot spot in the radio maps) to remain well defined during the evolution. Their morpholo- 



gies resemble those observed in naked quasar jets like 3C273 |36 ]. Fig. |T| shows several 
snapshots of the time evolution of a light, relativistic jet with large internal energy. The 
dependence of the beam's internal structure on the flow speed suggests that relativistic 
effects may be relevant for the understanding of the difference between slower, knotty BL 



Lac jets and faster, smoother quasar jets |60] 



Highly supersonic models, in which kinematic relativistic effects due to high beam 
Lorentz factors dominate, have extended overpressured cocoons. These overpressured co- 
coons can help to confine the jets during the early stages of their evolution | |107| and even 
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Figure 12: Time evolution of a light, relativistic (beam flow velocity equal to 0.99) jet with 
large internal energy. The logarithm of the proper rest-mass density is plotted in grey scale the 
maximum value corresponding to white and the minimum to black. 
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Figure 13: Logarithm of the proper rest-mass density and energy density (from top to bottom) 
of an evolved, powerful jet propagating through the intergalactic medium. The white contour 
encompasses the jet material responsible for the synchrotron emission. 
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Figure 14: Computed radio maps of a compact relativistic jet showing the evolution of a su- 
per luminal component (from left to right). Two resolutions are shown: present VLBI resolution 
(white contours) and resolution provided by the simulation (black/white images). 



cause their deflection when propagating through non-homogeneous environments ||145 
The cocoon overpressure causes the formation of a series of oblique shocks within the beam 
in which the synchrotron emission is enhanced. In long term simulations (see Fig. |13|), the 
evolution is dominated by a strong deceleration phase during which large lobes of jet ma- 
terial (like the ones observed in many FRIIs, e.g., Cyg A |24|| ) start to innate around the 
jet's head. These simulations reproduce some properties observed in powerful extragalactic 
radio jets (lobe inflation, hot spot advance speeds and pressures, deceleration of the beam 
flow along the jet) and can help to constrain the values of basic parameters (such as the 
particle density and the flow speed) in the jets of real sources. 

The development of multidimensional relativistic hydrodynamic codes has allowed, for 
the first time, the simulation of parsec scale jets and superluminal radio components 



84, 115 



The presence of emitting flows at almost the speed of light enhances the 
importance of relativistic effects in the appearance of these sources (relativistic Doppler 
boosting, light aberration, time delays). Hence, one should use models which combine 
hydrodynamics and synchrotron radiation transfer when comparing with observations. In 
these models, moving radio components are obtained from perturbations in steady rela- 
tivistic jets. Where pressure mismatches exist between the jet and the surrounding at- 
mosphere reconfinement shocks are produced. The energy density enhancement produced 
downstream from these shocks can give rise to stationary radio knots as observed in many 
VLBI sources. Superluminal components are produced by triggering small perturbations 
in these steady jets which propagate at almost the jet flow speed. One example of this 
is shown in Fig.[l4| (see also f68|), where a superluminal component (apparent speed 
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times the speed of light) is produced from a small variation of the beam flow Lorentz factor 
at the jet inlet. The dynamic interaction between the induced traveling shocks and the 
underlying steady jet can account for the complex behavior observed in many sources |p7 |. 



The first magnetohydrodynamic simulations of relativistic jets have been already un- 
dertaken in 2D |81], |8(| and 3D ||125| , |126|| to study the implications of ambient magnetic 



fields in the morphology and bending properties of relativistic jets. However, despite the 
impact of these results in specific problems like ,e.g., , the understanding of the misaligne- 
ment of jets between pc and kpc scales, these 3D simulations have not addressed the effects 
on the jet structure and dynamics of the third spatial degree of freedom. This has been 
the aim of the work undertaken by Aloy et al. . 

Finally, Koide et al. |[!| have developed a general relativistic MHD code and applied 



it to the problem of jet formation from black hole accretion disks. Jets are formed with a 
two-layered shell structure consisting of a fast gas pressure driven jet (Lorentz factor m 2) 
in the inner part and a slow magnetically driven outflow in the outer part both of which 
being collimated by the global poloidal magnetic field penetrating the disk. 

7.2 Gamma- Ray Bursts (GRBs) 

A second phenomenon which involves flows with velocities very close to the speed of light 
are gamma-ray bursts (GRBs). Although known observationally since over 30 years, until 
recently their distance ("local" or "cosmological" ) has been, and their nature still is a 
matter of controversial debate [[57], |112| , 140 , 141 . GRBs do not repeat except for a few 



soft gamma-ray repeaters. They are detected with a rate of about one event per day, and 
their duration varies from milliseconds to minutes. The duration of the shorter bursts and 
the temporal substructure of the longer bursts implies a geometrically small source (less 
than ~ c • 1msec ~ 100 km), which in turn points towards compact objects, like neutron 
stars or black holes. The emitted gamma-rays have energies in the range 30keV to 2 MeV. 
Concerning the distance of GRB sources major progress has occured through the obser- 
vations by the BATSE detector on board the Compton Gamma-Ray Observatory (GRO), 
which have proven that GRBs are distributed isotropically over the sky |[L11|| . Even more 



important the detection and the rapid availability of accurate coordinates (~ arc minutes) 
of the fading X-ray counterparts of GRBs by the BeppoSAX spacecraft beginning in 1997 



f33| , |143|| has allowed for subsequent successful ground based observations of faint GRB 
afterglows at optical and radio wavelength. In case of GRB 990123 the optical, X-ray and 
gamma-ray emission was detected for the first time almost simultaneously (optical obser- 
vations began 22 seconds after the onset of the GRB) |21], |]]. From optical spectra thus 
obtained redshifts of several gamma-ray bursts have been determined, e.g., GRB 970508 
(z = 0.835 [|1|, P71D , GRB 971214 (z = 3.42 g]), GRB 980703 (z = 0.966 g(J) and 
GRB 990123 (1.60 < z < 2.05 §), which confirm that (at least some) GRBs occur at 
cosmological distances. Assuming isotropic emission the inferred total energy of cosmolog- 
ical GRBs emitted in form of gamma-rays ranges from several 10 51 erg to 3 10 53 erg (for 
GRB 971214) gj|, and exceeds 10 54 erg for GRB 990123 g, ||. Updated information on 
GRBs localized with BeppoSAX, BATSE/RXTE(PCA) or BATSE/RXTE(ASM) can be 
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obtained from a web site maintained by Greiner |7T[ . 

The compact nature of the GRB source, the observed flux and the cosmological dis- 
tance taken together imply a large photon density. Such a source has a large optical depth 
for pair production. This is, however, inconsistent with the optically thin source indicated 
by the non-thermal gamma-ray spectrum, which extends well beyond the pair produc- 
tion threshold at 500 keV. This problem can be resolved by assuming an ultra-relativistic 
expansion of the emitting region, which eliminates the compactness constraint. The bulk 
Lorentz factor required are then W > 100 (see, e.g., ||141]] ) 

In April 1998 the pure cosmological origin of GRBs was challenged by the detec- 
tion of the Type Ib/c supernova SN 1998bw |TT], |62|] within the 8 arc minute error box 



of GRB 980425 |161| , |138|| . Its explosion time is consistent with that of the GRB, and 
relativistic expansion velocities are derived from radio observations of SN1998bw W7\ . 
BeppoSAX detected two fading X-ray sources within the error box, one being position- 
ally consistent with the supernova and a fainter one not consistent with the position of 
SN1998bw [ |138|| . Taken together these facts suggest a relationship between GRBs and 
SNelb/c, i.e., core collapse supernovae of massive stellar progenitors which have lost their 
hydrogen and helium envelopes pi UjL M. As the host galaxy ESQ 184-82 of SN 1998bw 



is only at a redshift of z — 0.0085 ||170|| and as GRB 980425 was not extraordinarily bright, 



GRB-supernovae are more than four orders of magnitude fainter (Etot^ = 7 10 47 erg for 



GRB 980425 |25|| ) than a typical cosmological GRB. However, the observation of the second 
fading X-ray source within the error box of GRB 980425 and unrelated with SN 1998bw 
still causes some doubts on the GRB supernova connection, although the probability of 
chance coincidence of GRB 980425 and SN 1998bw is extremely low |188| . 

In order to explain the energies released in a GRB various catastrophic collapse events 
have been proposed including neutron-star /neutron-star mergers [ |131| , |B9] , pE8] |, neutron- 
star /black-hole mergers |117|| , collapsars [|1 87| , |98| and hypernovae |132 |. These models all 
rely on a common engine, namely a stellar mass black hole which accretes several solar 
masses of matter from a disk (formed during a merger or by a non-spherical collapse) 
at a rate of 
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A fraction of the gravitational binding energy released by 
accretion is converted into neutrino and anti-neutrino pairs, which in turn annihilate into 
electron-positron pairs. This creates a pair fireball, which will also include baryons present 
in the environment surrounding the black hole. Provided the baryon load of the fireball is 
not too large, the baryons are accelerated together with the e + e~ pairs to ultra-relativistic 
speeds with Lorentz factors > 10 2 |26|, 142 , 141] . The existence of such relativistic flows is 
supported by radio observations of GRB 980425 [|87 |. It has been further argued that the 
rapid temporal decay of several GRB afterglows is inconsistent with spherical (isotropic) 
blast wave models, and instead is more consistent with the evolution of a relativistic jet 
after it slows down and spreads laterally ||156f| . Independent of the flow pattern the bulk 
kinetic energy of the fireball then is thought to be converted into gamma-rays via cyclotron 



radiation and/or inverse Compton processes (see, e.g., ||11I 

One-dimensional numerical simulations of spherically symmetric relativistic fireballs 
have been performed by several authors to model GRB sources ||142|, |133|, |134|. Multi- 



dimensional modelling of ultra-relativistic jets in the context of GRBs has for the first 
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time been attempted by Aloy et al. 0]. Using a collapsar progenitor model of MacFadyen 
& Woosley |9£| they have simulated the propagation of an axisymmetric jet through the 



mantle and envelope of a collapsing massive star (10 M©) using the GENESIS special 
relativistic hydrodynamic code 0. The jet forms as a consequence of an assumed energy 
deposition of 10 51 erg/sec within a 30 degree cone around the rotation axis. At break-out, 
i.e., when the jets reaches the surface of the stellar progenitor, the maximum Lorentz factor 
of the jet flow is about 20. 

8 CONCLUSION 

8.1 Evaluation of the methods 

An assessment of the quality of the numerical methods should consider, at least, the 
following aspects: (i) accuracy and robustness in describing high Lorentz factor flows with 
strong shocks; (ii) effort required to extend to multi dimensions; (iii) effort required to 
extend to RMHD and GRHD. In Table flT0| ) we have summarized these aspects of numerical 
methods for SRHD. 

Since their introduction in numerical RHD at the beginning of nineties, HRSC meth- 
ods have demonstrated their ability to describe accurately (stable and without excessive 
smearing) relativistic flows of arbitrarily large Lorentz factors and strong discontinuities 
reaching the same quality as in classical hydrodynamics. In addition (as it is the case for 
classical flows, too), HRSC methods show the best performance compared to any other 
method (e.g., artificial viscosity, FCT or SPH). 

Despite of the latter fact, a lot of effort has been put into improving these non HRSC 
methods. Using a consistent formulation of artificial viscosity has significantly enhanced 
the capability of finite difference schemes [ |128|| as well as of relativistic SPH |160|| to handle 



strong shocks without spurious post-shock oscillations. However, this comes at the price of 
a large numerical dissipation at shocks. Concerning relativistic SPH recent investigations 
using a conservative formulation of the hydrodynamic equations [2~9, |16U| | have reached an 



unprecedented accuracy with respect to previous simulations, although some issues still re- 
main. Besides the strong smearing of shocks, the description of contact discontinuities and 
of thin structures moving at ultrarelativistic speeds needs to be improved (see Sect. |6.2|). 

Concerning FCT techniques, those codes based on a conservative formulation of the 
RHD equations have been able to handle relativistic flows with discontinuities at all flow 
speeds, although the quality of the results is below that of HRSC methods in all cases 

The extension to multi dimensions is simple for most relativistic codes. Finite dif- 
ference techniques are easily extended using directional splitting. Note, however, that 



HRSC methods based on exact solutions of the Riemann problem [106, 182] first require 



the development of a multidimensional version of the relativistic Riemann solver. The 



adapt ing-gr id, artificial viscosity, implicit code of Norman & Winkler |128|| and the rela- 



tivistic Glimm method of Wen et al. ||182j| are restricted to one dimensional flows. Note 
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that Glimm's method produces the best results in all the tests analyzed in Sect.[| 



The symmetric TVD scheme proposed by Davis [37| and extended to GRMHD (see 



below) by Koide et al. [31] combines several characteristics making it very attractive. It is 
written in conservation form and is TVD, i.e., it is converging to the physical solution. In 
addition, it is independent of spectral decompositions, which allows for a simple extension 
to RMHD. Quite similar statements can be made about the approach proposed by van 
Putten ||175|| . In contrast to FCT schemes (which are also easily extended to general 



systems of equations), both Koide et al. 's and van Putten's method are very stable when 
simulating mildly relativistic flows (maximum Lorentz factors ~ 4) with discontinuities. 
Their only drawback is an excessive smearing of the latter. A comparison of Davis' method 
with Riemann solver based methods would be desirable. 

8.2 Further developments 

The directions of future developments in this field of research are quite obvious. They can 
be divided into four main categories: 

8.2.1 Incorporation of realistic microphysics 

Up to now most astrophysical SRHD simulations have assumed matter whose thermody- 
namic properties can be described by an inviscid ideal equation of state with a constant 
adiabatic index. This simplification may have been appropriate in the first generation 
of SRHD simulations, but it clearly must be given up when aiming at a more realistic 
modelling of astrophysical jets, gamma-ray burst sources or accretion flows onto compact 
objects. For these phenomena a realistic equation of state should include contributions 
from radiation (7 = 4/3 "fluid"), allow for the formation of electron-positron pairs at 
high temperatures, allow the ideal gas contributions to be arbitrarily degenerate and/or 
relativistic. Effects due to cooling, heat conduction, nuclear transmutations and viscosity 
may have to be considered, too. When simulating relativistic heavy ion collisions the use 
of a realistic equation of state is essential for an adequate description of the phenomenon. 
However, as these simulations have been performed with FCT based difference schemes 
(see, e.g., |162| ), this poses no specific numerical problem. The simulation of flows obeying 



an elaborated microphysics with HRSC methods needs in some cases the extension of the 
present relativistic Riemann solvers to handle general equations of state. This is the case 
of the Roe-Eulderink method (extendable by the procedure developed in the classical case 
by Glaister |64]]), and rPPM and rGlimm both relying on a exact solution of the Riemann 
problem for ideal gases with constant adiabatic exponent (which can also be extended fol- 
lowing the procedure of Colella & Glaz |JT| for classical hydrodynamics). We expect the 



second generation of SRHD codes being capable of treating general equations of state and 
various source/sink terms routinely. 



52 



Table 10: Evaluation of numerical methods for SRHD. Methods have been categorized for clarity 



Extension to several 
spatial dimensions b 



Method 



Ultrarelativistic 
regime 



Handling of 
discontinuities a 



Extension to 
GRHD RMHD 



AV-mono 


X 


0, SE 


V 


V 


X 


cAV-implicit 


V 


V 


X 


X 


X 


HRSC C 


V 


V 


v d 


v e 


x / 


rGlimm 


V 


V 


X 


X 


X 


sTVD 


V 9 


D 


V 


v 7 


V 


van Putten 


V 9 


D 


V 


X 


V 


FCT 


V 





V 


X 


X 


SPH 


V 


D, 


V 


v h 


X l 



a D: excessive dissipation; O: oscillations; SE: systematic errors. 

b all finite difference methods are extended by directional splitting. 

c contains all the methods listed in Table with exception of rGlimm [[182|| and sTVD HSl| . 

d rPPM [p.06|] requires an exact relativistic Riemann solver with non-zero transverse 



GRHD extensions of several HRSC methods based on linearized Riemann 



speeds. 

e there exist 

solvers. The procedure developed by Pons et al. ||146| allows any SRHD Riemann solver to 

be applied to GRHD flows. 

* except HLL which requires spectral decomposition of RMHD equations or solution of 
RMHD Riemann problem. Van Putten |180|| has studied the characteristic structure of 
the RMHD equations in (constraint free) divergence form as a first step to extend modern 
HRSC methods to RMHD. Komissarov [Q has developed a multidimensional RMHD 
code based on a linearized Riemann solver. 

9 needs confirmation. 
h codes of refs. [|^, |160| . 

* There is one code which considered such an extension |J9], but the results are not 
completely satisfactory. 
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8.2.2 Coupling of SRHD schemes with AMR 

Modelling astrophysical phenomena often involves an enormous range of length scales and 
time scales to be covered in the simulations. In two and definitely in three spatial dimen- 
sions many such simulations cannot be performed with sufficient spatial resolution on a 
static equidistant or non-equidistant computational grid, but they will require dynamic, 
adaptive grids. In our opinion the most promising approach in this direction will be the 
coupling of SRHD solvers with the adaptive mesh refinement (AMR) techniques [|12 |. AMR 
automatically increases the grid resolution near flow discontinuities or in regions of large 
gradients (of the flow variables) until a prescribed accuracy of the difference approxima- 
tion is achieved. A SRHD simulation of a relativistic jet based on a combined HLL-AMR 
scheme was performed by Duncan & Hughes EIJ . Plewa et al. ||145| have modelled the 



deflection of highly supersonic jets propagating through non-homogeneous environments 



AMRA of Plewa 

with the adaptive grid code 

for industrial applications 



using the HRSC scheme of Marti et al. [ |108|| combined with the AMR implementation 
144 ]. Komissarov & Falle [|84| have combined their numerical scheme 
Cobra, which has been developed by Mantis Numerics Ltd. 
and which uses a hierarchy of grids with a constant refine- 



ment factor of two between subsequent grid levels. 

8.2.3 General relativistic hydrodynamics (GRHD) 

Up to now only very few attempts have been made to extend HRSC methods to GRHD 
and all of these have used linearized Riemann solvers ||104| , 51 , |153| , H |59|. In the most 
recent of these approaches Font et al. 
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have developed a 3D general relativistic HRSC 
hydrodynamic code where the matter equations are integrated in conservation form and 
fluxes are calculated with Marquina's formula. A very interesting and powerful procedure 
has been proposed by Pons et al. [|146|1 , which allows one to exploit all the developments 
in the field of special relativistic Riemann solvers in general relativistic hydrodynamics. 
The procedure relies on a local change of coordinates at each zone interface such that 
the spacetime metric is locally flat. In that locally flat spacetime any special relativistic 
Riemann solver can be used to calculate the numerical fluxes, which are then transformed 
back. Finer grids and improved time advancing methods will be required in regions where 
large gradients or temporal variations of the gravitational field are encountered. The 
numerical implementation is simple and computationally inexpensive. 

The characteristic formulations of the Einstein field equations are able to handle the 
long term numerical description of single black hole space times in vacuum |TJ]]. In order 
to include matter in such an scenario, Papadopoulos & Font ||135|1 have generalized HRSC 
methods to cope with the hydrodynamic equations in this null foliation. 

Other developments in GRHD in the past included finite element methods for simulat- 
ing spherically symmetric collapse in general relativity ||101|| , general relativistic pseudo- 
spectral codes based on the (3+1) ADM formalism Q for computing radial perturbations 
fTOfl and 3D gravitational collapse of neutron stars |Tj| , and general relativistic SPH |)9| . 



The potential of these methods for the future is unclear, as none of them is specifically 
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appropriate for ultrarelativistic speeds and strong shock waves which are characteristic of 
most astrophysical applications. 
Peitz & Appl 



136 



have addressed the difficult issue of non-ideal GRHD, which is of 
particular importance, e.g., for the simulation of accretion discs around compact objects, 
rotating relativistic fluid configurations, and the evolution of density fluctuations in the 
early universe. They have accounted for dissipative effects by applying the theory of 
extended causal thermodynamics, which eliminates the causality violating infinite signal 
speeds arising from the conventional Navier-Stokes equation. Peitz & Appl have not 
implemented their model numerically yet. 



8.2.4 Relativistic magneto— hydrodynamics (RMHD) 

The inclusion of magnetic effects is of great importance in many astrophysical flows. The 
formation and collimation process of (relativistic) jets most likely involves dynamically 
important magnetic fields and occurs in strong gravitational fields. The same is likely to 
be true for accretion discs around black holes. Magneto-relativistic effects even play a non 
negligible role in the formation of proto-stellar jets in regions close to the light cylinder 
22[ . Thus, relativistic MHD codes are a very desirable tool in astrophysics. The non-trivial 



task of developing such a kind of code is considerably simplified by the fact that because 
of the high conductivity of astrophysical plasmas one must only consider ideal RMHD in 
most applications. 

Special relativistic 2D MHD test problems with Lorentz factors up to ~ 3 have been 
investigated by Dubai J44j with a code based on FCT tecniques (see Sect. ^). In a series of 

p5jp6|,[8^ 



papers Koide and coworkers J5I], ^, |125| , |126| , [82]] have investigated relativistic magnetized 
jets using a symmetric TVD scheme (see Sect. 0). Koide, Nishikawa & Mutel |81| simulated 



a 2D RMHD slab jet, whereas Koide |3(J investigated the effect of an oblique magnetic 
field on the propagation of a relativistic slab jet. Nishikawa et al. |125| , |126| extended 
these simulations to 3D and considered the propagation of a relativistic jet with a Lorentz 
factor W = 4.56 along an aligned and an oblique external magnetic field. The 2D and 3D 
simulations published up to now only cover the very early propagation of the jet (up to 20 
jet radii) and are performed with moderate spatial resolution on an equidistant Cartesian 
grid (up to 101 zones per dimension, i.e., 5 zones per beam radius). 

Van Putten [ |174j , |175| | has proposed a method for accurate and stable numerical simula- 
tions of RMHD in the presence of dynamically significant magnetic fields in two dimensions 
and up to moderate Lorentz factors. The method is based on MHD in divergence form 
using a 2D shock-capturing method in terms of a pseudo-spectral smoothing operator (see 
Sect. |. He applied this method to 2D blast waves [ 177 1 and astrophysical jets [ 176 , 179 1. 
Steps towards the extension of linearized Riemann solvers to ideal RMHD have already 
been taken. Romero |154|| has derived an analytical expression for the spectral decomposi- 
tion of the Jacobian in the case of a planar relativistic flow field permeated by a transversal 
magnetic field (nonzero field component only orthogonal to flow direction). Van Putten 
180| has studied the characteristic structure of the RMHD equations in (constraint free) 
divergence form. Finally, Komissarov p3) has presented a robust Godunov-type scheme for 
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RMHD, which is based on a linear Riemann solver, has second-order accuracy in smooth 
regions, enforces magnetic flux conservation, and which can cope with ultrarelativistic 
flows. 



We end with the simulations performed by Koide, Shibata & Kudoh [82| on magnet- 
ically driven axisymmetric jets from black hole accretion disks. Their GRMHD code is 
an extension of the special relativistic MHD code developed by Koide et al. |81|, [8(], |125| . 
The necessary modifications of the code were quite simple, because in the (nonrotating) 
black hole's Schwarzschild spacetime the GRMHD equations are identical to the SRMHD 
equations in general coordinates, except for the gravitational force terms and the geometric 
factors of the lapse function. With the pioneering work of Koide, Shibata & Kudoh the 
epoch of exciting GRMHD simulations has just begun. 

9 ADDITIONAL INFORMATION 

9.1 Algorithms to recover primitive quantities 

The expressions relating the primitive variables (p, v \ p) to the conserved quantities (D, S\ r) 
depend explicitly on the equation of state p(p, e) and simple expressions are only obtained 
for simple equations of state (i.e., ideal gas). 

A function of pressure, whose zero represents the pressure in the physical state, can 
easily be obtained from Eqs. (@-g3)> (0) and (0): 



f{p)=p{p*{p),£*(p))-P (61) 

with p*(p) and £*(p) given by 

P.(P)~, (62) 

and 

£ * {P) ~ DWM ' {b6) 

where 

WM = l (64) 

\/l-vl(p)v #i (p) 

and 

vi(p) = ^ • (65) 

* v ' T + D+p y ' 

The root of flBTD can be obtained by means of a nonlinear root-finder (e.g., a one-dimensional 
Newton-Raphson iteration). For an ideal gases with a constant adiabatic exponent such 
a procedure has proven to be very successful in a large number of tests and applications 
|104| , |106| , |108|| . The derivative of / with respect to p, /', can be approximated by 



f' = vl(p)v^p)c s ^p) 2 -l, (66) 
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where c s * is the sound speed which can efficiently be computed for any EOS. Moreover, 
approximation (^) tends to the exact derivative as the solution is approached. 

Eulderink [^, |51| has also developed several procedures to calculate the primitive 
variables for an ideal EOS with a constant adiabatic index. One procedure is based on 
finding the physically admissible root of a fourth-order polynomial of a function of the 
specific enthalpy. This quartic equation can be solved analytically by the exact algebraic 
quartic root formula although this computation is rather expensive. The root of the quartic 
can be found much more efficiently using a one-dimensional Newton-Raphson iteration. 
Another procedure is based on the use of a six-dimensional Newton-Kantorovich method 
to solve the whole nonlinear set of equations. 

Also for ideal gases with constant gamma, Schneider et al. [|157|| transform system (||- 



|T0D , ( |r2| ) and ( jl3|) algebraically into a fourth-order polynomial in the modulus of the flow 
speed, which can be solved analytically or by means of iterative procedures. 

For a general EOS, Dean et al. JJUJ and Dolezal & Wong |5]| proposed the use of iterative 
algorithms for v 2 and p, respectively. 

9.2 Spectral decomposition of the 3D SRHD equations 

The full spectral decomposition (right and left eigenvectors) of the Jacobian matrices as- 
sociated to the SRHD system in 3D has been first derived by Donat et al. p2| . Previously, 
Marti et al. . |104j| obtained the spectral decomposition in ID SRHD and Eulderink pOf 



and Font et al. [58|, the (eigenvalues and) right eigenvectors in 3D. The Jacobians are given 
by 

* = ^, (67) 

where the state vector u and the flux vector F* are defined in (||) and (|7|), respectively. 
In the following we explicitly give both the eigenvalues and the right and left eigenvec- 
tors of the Jacobi matrix B x only (the cases i = y,z are easily obtained by symmetry 
considerations) . 
The eigenvalues of matrix B x (u) are 

A± = - _ l ^ c2 |w x (l - c 2 s )±c s ^J(l - v 2 )[l - W* - {v 2 - v x v x )c 2 }} , (68) 

and 

\ = v x (triple). (69) 

A complete set of right-eigenvectors is 

r 02 = (Wv\ 2hW 2 v x v y , h{l + 2W 2 v y v y ), 2hW 2 v y v z , 2hW 2 v y - Wv y ) (71) 
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r 0i3 = (Wv z , 2hW 2 v x v z , 2hW 2 v y v z , h{\ + 2W 2 v z v z ), 2hW 2 v z - Wv' 



r± = (1, hWA±X±, hWv y , hWv z , hWA± - I) 



where 



K K . 1 — V X V x 

K ~ ~ 2 ' K = ~ ' A ±= i ST~ 

K — Cg p 1 — U X A± 

The corresponding complete set of left-eigenvectors is 

W 



lo,i 



/C-l 



(/i - W, Wu x , WV, WV, -W) 



(72) 
(73) 

(74) 



H),2 



/l(l — U^V 1 ) 



(-v y ,v x v y ,l-v x v x ,0,-v y ) 



lo,3 



-(—■ y z , f/V, 0, 1 — f V, — f 2 ) 



1 



/? 2 

:±1 4 



/l(l — v x l; x ) 
^^(v 1 - A±) -v x - W 2 (v 2 - v x v x ){2K - l)(v x - A±\±) + ICA ± X ± 

1 + W 2 {v 2 - v x v x )(2K - 1)(1 - At) - £.A± 

WV(2/C - l)A±(v x - A±) 

WV(2/C - l)At(u a - A±) 

-^ - W^ 2 (w 2 - uV)(2/C - l)(u x - A±\±) + /C4±A± 
where A is the determinant of the matrix of right-eigenvectors, i.e., 

A = h 3 W(}C - 1) (1 - v x v x ) (A+X+ - A-X.) . (75) 

For an ideal gas equation of state K, — h, i.e.,JC > 1, and hence A^O for \v x \ < 1. 

9.3 Program RIEMANN 

PROGRAM RIEMANN 

C C This program computes the solution of a ID 

c relativistic Riemann problem with 

C initial data UL if X<0.5 and UR if X>0.5 

C in the whole spatial domain [0, 1] 

END 
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9.4 Basics of HRSC methods and recent developments 

In this section we introduce the basic notation of finite differencing and summarize recent 
advances in the development of HRSC methods for hyperbolic systems of conservation 
laws. The content of this section is not specific to SRHD, but applies to hydrodynamics 
in general. 

In order to simplify the notation and taking into account that most powerful results 
have been derived for scalar conservation laws in one spatial dimension, we will restrict 
ourselves to the initial value problem given by the equation 

du df(u) 

m + ^ = ° (76) 

with the initial condition u(x, t = 0) = Uq(x). 

In hydrodynamic codes based on finite difference or finite volume techniques, equation 
( f76f) is solved on a discrete numerical grid (xj,t n ) with 

Xj = (j - 1/2) Ax, j = 1,2,..., (77) 

and 

t n = nAt, n = 0,1,2,... , (78) 

where At and Ax are the time step and the zone size, respectively. A difference scheme 
is a time-marching procedure allowing one to obtain approximations to the solution at 
the new time, w™ +1 , from the approximations in previous time steps. Quantity w™ is an 
approximation to u(xj,t n ) but, in the case of a conservation law, it is often preferable to 
view it as an approximation to the average of u(x, t) within a zone [3^-1/2,2^+1/2] (i.e., as 
a zone average), where Xj±i/2 — ( x j + x j±i)/2- Hence 

1 r x j+i/2 
u"= — / u{x,t n )dx, (79) 



Ax JXj-1/3 

which is consistent with the integral form of the conservation law. 

Convergence under grid refinement implies that the global error 1 1 E& x 1 1 , defined as 

\\E Ax \\=AxJ2\u]-u]\, (80) 

3 

tends to zero as Ax — > 0. For hyperbolic systems of conservation laws methods in con- 
servation form are preferred as they guarantee that if the numerical solution converges, it 
converges to a weak solution of the original system of equations (Lax-Wendroff theorem 
]). Conservation form means that the algorithm can be written as 



W„- " — Uj \f{Uj_ r ,Uj_ r+1 ,...,Uj +q ) j{Uj_ r _ 11 Uj_ r ,...,Uj +q _i)J 



where q and r are positive integers, and / is a consistent (i.e.,f(u,u,...,u) = /(«)) 
numerical flux function. 
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The Lax-Wendroff theorem cited above does not establish whether the method con- 
verges. To guarantee convergence, some form of stability is required, as for linear problems 
(Lax equivalence theorem |150|| ). In this context the notion of total-variation stability has 
proven to be very successful, although powerful results have only been obtained for scalar 
conservation laws. The total variation of a solution at t = t n , TV(u n ), is defined as 

+00 
TV(«") = £ \u] +1 - u]\ . (82) 

A numerical scheme is said to be TV-stable, if TV(u n ) is bounded for all At at any time 
for each initial data. One can then prove the following convergence theorem for non-linear, 
scalar conservation laws |95| : For numerical schemes in conservation form with consistent 
numerical flux functions, TV-stability is a sufficient condition for convergence. 

Modern research has focussed on the development of high-order, accurate methods in 
conservation form, which satisfy the condition of TV-stability The conservation form is 
ensured by starting with the integral version of the partial differential equations in conser- 
vation form (finite volume methods). Integrating the PDE over a finite space-time domain 
[ x j-i/2,Xj+i/2] x [t n ,t n+1 ] and comparing with fl8l|) , one recognizes that the numerical flux 
function fj+1/2 is an approximation to the time-averaged flux across the interface, i.e., 

1 r tn+1 
/i+i/2~^y <n f(u{x j+1/2 ,t))dt. (83) 

Note that the flux integral depends on the solution at the zone interface, u(xj + i/2, t), during 
the time step. Hence, a possible procedure is to calculate u(xj+i/2, t) by solving Riemann 
problems at every zone interface to obtain 

u(x j+1/2 , t) = w(0; u], u] +1 ) . (84) 

This is the approach followed by an important subset of shock-capturing methods, called 
Godunov-type methods J73], |U| after the seminal work of Godunov [BB], who first used 



an exact Riemann solver in a numerical code. These methods are written in conserva- 
tion form and use different procedures (Riemann solvers) to compute approximations to 
ti(0; u^ , Uj +l ) . The book of Toro ||171|| gives a comprehensive overview of numerical meth- 



ods based on Riemann solvers. The numerical dissipation required to stabilize an algorithm 
across discontinuities can also be provided by adding local conservative dissipation terms to 
standard finite-difference methods. This is the approach followed in the symmetric TVD 
schemes developed in |37|, |152| , |192|| . 



High-order of accuracy is usually achieved by using conservative monotonic polynomial 
functions to interpolate the approximate solution within zones. The idea is to produce more 
accurate left and right states for the Riemann problem by substituting the mean values 
u n - (that give only first-order accuracy) by better representations of the true flow near 
the interfaces, let say u^ +1 i 2) u^ +1 , 2 . The FCT algorithm [19| constitutes an alternative 
procedure where higher accuracy is obtained by adding an anti-diffusive flux term to 
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the first-order numerical flux. The interpolation algorithms have to preserve the TV- 
stability of the scheme. This is usually achieved by using monotonic functions which lead 
to the decrease of the total variation (total-variation-diminishing schemes, TVD |72[ ) . 
High-order TVD schemes were first constructed by van Leer ||172 |, who obtained second- 
order accuracy by using monotonic piecewise linear slopes for cell reconstruction. The 
piecewise parabolic method (PPM) J32] provides even higher accuracy. The TVD property 
implies TV-stability, but can be too restrictive. In fact, TVD methods degenerate to 
first-order accuracy at extreme points ||130|. Hence, other reconstruction alternatives have 



been developed where some growth of the total variation is allowed. This is the case for 
the total-variation-bounded (TVB) schemes ||158| , the essentially non-oscillatory (ENO) 
schemes ]73| and the piecewise-hyperbolic method (PHM) ||102| . 



9.5 Newtonian SPH equations 



Following Monaghan ||120| the SPH equation of motion for a particle a with mass m and 
velocity v is given by 



dVa 

dt 



-J2 m >> 



A pi J 



*5) 



where the summation is over all particles other than particle a, p is the pressure, p is the 
density, and d/dt denotes the Lagrangian time derivative. II a & is the artificial viscosity 
tensor, which is required in SPH to handle shock waves. It poses a major obstacle in 
extending SPH to relativistic flows (see, e.g., |77|, |29|). W ab is the interpolating kernel, and 
V ' a W ab denotes the gradient of the kernel taken with respect to the coordinates of particle 
a. 

The kernel is a function of |r — r&| (and of the SPH smoothing length /isph)> i.e., its 
gradient is given by 

V a W a6 = v ab F ab , (86) 

where F ab is a scalar function which is symmetric in a and b, and r ab is a shorthand for 
(r a — 17,). Hence, the forces between particles are along the line of centers. 

Various types of spherically symmetric kernels have been suggested over the years |118 



11 1 . Among those the spline kernel of Monaghan & Lattanzio ||121|| , mostly used in current 
SPH-codes, yields the best results. It reproduces constant densities exactly in ID, if the 
particles are placed on a regular grid of spacing h S p H , and has compact support. 
In the Newtonian case II afe is given by [|120| 



n 



ab 



hsPH v a fe • r ab I _ 

-« — — i T5 Cab 

Pab VabV V 



^sph v a & • r ab 

\ r ab\ 2 



m 



provided v ab ■ r ab < 0, and II ab = otherwise. Here v a6 = v a - v b , c ab = ~(c a + c b ) is the 
average sound speed, p ab = \{p a + Pb), and a ~ 1.0 is a parameter. 
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Using the first law of thermodynamics and applying the SPH formalism one can derive 



the thermal energy equation in terms of the specific internal energy e (see, e.g., ||119f[ ). How- 
ever, when deriving dissipative terms for SPH guided by the terms arising from Riemann 
solutions, there are advantages to use an equation for the total specific energy E = v 2 /2+e, 



which reads [120] 



dE a „ (p a \ b p b v a 



^ m b ^ + ^ + il ab ■ V„W„ 6 , (88) 



where Q ab is the artificial energy dissipation term derived by Monaghan [120|. For the 
relativistic case the explicit form of this term is given in Section fl4.2|) . 

In SPH calculations the density is usually obtained by summing up the individual 
particle masses, but a continuity equation may be solved instead, which is given by 

^7 = - E Mv a - v b ) V a W ab . (89) 

at b 



The capabilities and limits of SPH have been explored, e.g., in |164| , |167|| . Steinmetz 



& Miiller ||164| | conclude that it is possible to handle even difficult hydrodynamic test 



problems involving interacting strong shocks with SPH provided a sufficiently large number 
of particles is used in the simulations. SPH and finite volume methods are complementary 
methods to solve the hydrodynamic equations each having its own merits and defects. 
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